Spectral Learning of Shared Dynamics Between Generalized-Linear Processes
Authors: Lucine L Oganesian
Omid G. Sani
Maryam Shanechi
Venue: NeurIPS 2024 poster (2024)
Keywords: state space models
subspace identification
dynamical systems
neural coding
generalized-linear models
Primary area: neuroscience_and_cognitive_science

Abstract: Generalized-linear dynamical models (GLDMs) remain a widely-used framework within neuroscience for modeling time-series data, such as neural spiking activity or categorical decision outcomes. Whereas the standard usage of GLDMs is to model a single data source, certain applications require jointly modeling two generalized-linear time-series sources while also dissociating their shared and private dynamics. Most existing GLDM variants and their associated learning algorithms do not support this capability. Here we address this challenge by developing a multi-step analytical subspace identification algorithm for learning a GLDM that explicitly models shared vs. private dynamics within two generalized-linear time-series. In simulations, we demonstrate our algorithm's ability to dissociate and model the dynamics within two time-series sources while being agnostic to their respective observation distributions. In neural data, we consider two specific applications of our algorithm for modeling discrete population spiking activity with respect to a secondary time-series. In both synthetic and real data, GLDMs learned with our algorithm more accurately decoded one time-series from the other using lower-dimensional latent states, as compared to models identified using existing GLDM learning algorithms.

Spectral Learning of Shared Dynamics Between Generalized-Linear Processes
Lucine L. Oganesian
Ming Hsich Department of Electrical and Computer Engineering University Southern California
Los Angeles, CA
loganesiQusc.edt
Omid G. Sani
Ming Hsieh Department of Electrical and Computer Engineering
University of Southern California
Los Angeles, CA
omid. gwwsemsaaniiiiiio antmusc..iuc....m edu
Maryam M. Shanechi
Ming Hsieh Department of Electrical and Computer Engineering
Thomas Lord Department of Computer Science
Alfred E. Mann Department of Biomedical Engineering Neuroscience Graduate Program University Southern California
Los Angeles, CA shanech @usc. edu
Abstract
Generalized-linede dynamical models (GLDMs) remain a widely-used framework within neuroscience for modeling time-series data, such neural spiking activity categorical decision outcomes. Whereas the standard usage of GLDMs is to model single data source, certain applications require jointly modeling linear time-series sources while also dissociating their shared and private dynamics. Most existing GLDM variants and their associated learning algorithms do not support this capability. Here we address this challenge by developing multi-step analytical subspace identification algorithm for learning a GLDM that explicitly models shared vs. private dynamics within two generalized-linee time-scries. In simulations, we demonstrate our algorithm's s ability to dissociate and model the dynamics within two time-series sources while being agnostic to their respective observation distributions. In neural data, we consider two specific applications of our algorithm for modeling discrete population spiking activity with respect to a secondary time-series. both synthetic and real data, GLDMs learned with our algorithm more accurately decoded one time-series from the other using lower- dimensional latent states, as compared models identified using existing GLDM learning algorithms.
generalized-
"Corresponding author. Project code: http://ghhab..oomSShhchhhhlllllu
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Introduction
Generalized- linear dynamical models (GLDMs) are commonly used framework for modeling dynamics using a -ow-dimensional latent variable that evolves over time 1133. Due to their inter- pretability, data efficiency, and amenability to real-time engineering operations, GLDMs remain widely popular tool neuroscience for modeling time-series data, whether it be Poisson spiking neural activity. Bernoulli/Bnomiial categorical task variables, Gaussian behavior 49. Whereas most existing GLDM variants and their associated analytical learning algorithms focus on modeling the dynamics within a single time-series, there exist applications that require explicit dissociation of shared vs. private dynamics within two generalized-linen observation sources. For example, such functionality is helpful when modeling the dynamical relationship between recorded neural activity and certain behaviors of interest [10H55]
Here we fill these methodological gaps by deriving a novel covariance- baated subspace system identification (SSID) algorithm that is capable, with its multi-staged learning approach, of identifying the shared dynamics between two eeneralized-linear time-series with priority, before modeling the dynamics private to each observation. We design the method to seamlessly generalize to different observation distributions, such as Poisson Bernoulli. To illustrate the method, first show in simulations that our method successfully dissociates the shared dynamics within two generalized- linear time-series, agnostic of their respective observation models; compare against existing GLDM methods, we focused on Poisson, Bernoulli, and Gaussian generalized- observations. Next we demonstrate our method on two public non-human primate (NHP) datasets of discrete population spiking activity recorded from different brain regions and during different contexts [ 18. Compared with existing Poisson GLDMs and their learning algorithms, our method learned models that more accurately decoded one time-series from the other using lower-dimensional latent states, suggesting improved learning of shared dynamics.
2 Background
For linear state-space models with continuous Gaussian observations, subspace system identification (SSID) theory provides computationally efficient non-iterative algorithms for analytically learning state-space models, both with and without identification of shared dynamics 14 14]. These methods, however, either are not applicable to generalized-linett time-series with Gaussaan observations [14 19 20 23 or do not have the ability to dissociate shared vs. private dynamics between two time-series 21 22 24
To help with the exposition of our method in section [3] first review standard covariance- basee SSID and an existing SSID method for modeling Poisson point-processes 21, a widely-used class of generalized-linet observations.
2.1 Standard covariance-based SSID
The standard formulation for a linear state-space model with continuous Gaussian observations is as
{
X = Ax + We
(I)
rk
C,X² V
where X R is the latent state variable, rς R corresponds continuous Gaussian ob- servations, and w and v are state and observation noise terms, respectively, with distributions N(w;; 0, Q) and N(v 0, R), and cross-covariance S. Further, we define G == Cov (xk-x rk),, as the covariance between future latent state and current observation, and Ar. == Cov rk), as the instantaneous covariance the observations. Standard covariance- based SSID learns the parameters of a latent dynamical system θ (A, Cr. Ar.) given training samples т₁ and hyperparameter n₂ that specifies the latent state dimensionality. To do so, a future-paa past Hankel matrix, Hr, is first
constructed from the cross-covariances of the system's 's linear observations as [19 20
Ar:
Ar.-
Ar
Ariaa
Ar.
A
F;
Fo
H, Co)(
Tp
:
Arz-n
Ar
r2
where the integer denotes the user-specified maximum temporal lag (i.e., horizon) used to construct H, and Ar, := Cov rk) the 7-th lag cross-covariance for any timepoint k; under time- stationary assumptions. We note that the rank of Hr must be at least ng in order to identify

---

model with a latent dimension of ng. Thus the user-specified horizon must satisfy i X n, 2 n2g- Covariance- based SSID then decomposes H, into a product of observability (I,) and controllability
(A) matrices as 19 20
C₁ C,A
H, A
[A-
AG G]
(3)
C,A
where G is defined as above. The factorization of H, is done by computing a singular value decomposition (SVD) of Hr and keeping the top nhr singular values and corresponding singular vectors. From the factors of Hr, C, is read off as the first ny rows of T and G' is read off as the last ny columns of A. A is learned by solving r, A,A, where H, and D denote r, from which the top or bottom nh rows have been removed, respectively. This optimization problem has the following closed-form least-squares solution A = r.rr. with denoting the pseudo-inverse operation. The final parameter Aro is computed the empirical covariance of rg. See appendix A.4 on how (G, Ar.) specify (Q, R. S).
2.2 SSID for a single generalized -lnnear time-series
There has been some work extending SSID to generalized-lineat time-series, such as Poisson and Bernoulli observations 21 24 These methods, however, only learn the dynamics of single generalized- linear time-series rather than model shared vs. private dynamics between two time-series. Here we present one of our baselines, PLDSID [2I], which models a single Poisson series, as an example. A Poisson linear dynamical system (PLDS) model is defined as
X1+1 rk = C,X: + Yk rk
Ax + Wa
Poisson(exp(rr))
where X; E IR= is the latent state as before and Yı E RR= corresponds to discrete (e.g., neural spiking) observations which, conditioned on the latent process rı. are Poisson-distribueed with rate equal to the exponential rę (i.e., log-rate). Finally, W, Gtaussianddistributer state noise with covariance parameter Q, as before, and b is constant baseline log-rate. The PLDS model is commonly used for modeling Poisson process events, such as neural spiking activity 12 4 21 25
Buesing al. (21) developed SSID algorithm, termed PLDSID, to learn the PLDS model parameters OPLDS = (A, Cr: b, Q) given training samples yı and hyperparameter T2z-
SSID algorithms (section 2.1 are not directly applicable to Poisson- distributed observations. This is because the log-rates F& that are linearly related the latent states in equation 4 are not observable in practice rather, only a stochastic Poisson emission from them (i.e., Y&) is observed. As a result, the second moments constituting H, (i.e., Ar. ) cannot be directly estimated. The critical insight by Buesing et al. [21] was to leverage the log link function (i.e., exp 11) and the known conditional distribution yı F to compute the first (ft+=) ) and second (A,+) ) moments of the log-rate F3 from the first (My=1 ) and second (Ay=) moments of the discrete observations Yи. The + denotes that moments are computed for the future-past stacked vector of observations r [r] rp .T
Standard covariance-
T and y [y] yT] where
Hp+ : Err+] My+ == Ely+1 Apt == Cov(r-, r Ayt := Cov(y y=*).
To compute moments of the log-rate, Buesing al. (21] derived the following moment conversion
P mm
+ u12 - -My ym In(Ay mi.. + +2² y -My (In(42
2 11n
2 In(A, ymm
(5)
P
A + = In y + Myyt Myy) In(py My
where m A n correspond to different indices of the first and second moments of the future-past stacked observation vectors r and y=, and n, = =I.... Kn, where K is the total number of time points. With the first and second moments computed the moment conversion above, the baseline log rate b parameter is read off the first rows of μet and the Hankel matrix, H,, is constructed as per equation
(2). From here, we can proceed with the standard covariance-based SSID algorithm using Hr, as outlined in section 2.1 Discussion regarding learning the state noise covariance parameters (e.g., Q) is postponed to appendix section A.3) where we use an approach that unlike Buesing et al. 21 ensures validity of learned noise statistics.
3

---

Method
3.1 Model definition and assumptions
Both the standard linear state-space model (equation (T) and PLDS model (equation only model single observation on its own. To enable identification shared and private dynamics between two generalized-linee time-series, we write the following general multi-observation GLDM
Xk-1
Ax Wh
C,Su + Va +
Zk
C2x E₁
(6)
Yk
Pyir(ya: g(xk) Peiz(te; h(*)))
where Xg, rk: Wk. Vk. and b are defined in equations and to represent the second generalized- construction) to represent the latent process underlying this second observation, E ~ ~N(c; 0, F)) to represent the associated noise term, and d to represent the associated baseline value. We generically denote the probability distribution for уı conditioned PLDS model Pyir == (exp(on(ee((r))))) ) Pttz is defined similarly but for t. and 2g. Finally, g(:) and h(-) correspond to the link function the generalized- linear model, for example g(r;) = exp(rı) in PLDS models. In order dissociate between shared and private dynamics within the observation time-series, we introduce the following definition:
I We introduce observation time-series, 23 E R (with n4 = A2 by
the latent rı with Py. For example, in the
Definition 3.1. We take the system be written in block structure form defined below allowing us to dissociate shared from private latents
An
0
(1)
A
A2
A 22
0
C₂
= c (4)
c
C, = c
0
x
x
333
(3)
where X3
(1) Rn corresponds to latent states that drive both 2д and rg. X
(2)
I corresponds to
states that only drive rk. and (3) E TR's corresponds to states that only drive 2μ- with total states
nr n1 12 + ng. The parameter G can also be written in block partition format such that
E r E r
E E(ra]" IT
G =
X
-
ElraJ"
E (2)
E
G
E (3) F
E 3
E
G
can further simplify the definition of using following assumptions:
Assumption 3.2. The state noise covariance is assumed to have block diagonal structure such that Q = diag(Q) (1.2) Q(3)), where Q a square matrix of dimension nı + 122 and Q(3) is square matrix of dimension T23. Formally, the superseript notation (1,2) designates attribution the parameter to the first and second set of latent states.
Assumption 3.3. Initial latent states are assumed to be maually--ndeeendee making Cov Xn) diagonal.
These assumptions allow us to fully decouple the private latent states of the secondary time-series (x) from the latent states driving the primary time-series (x(2), (())). As result, can take (3) = (ov(x 1 r&) = 0. From the perspective state estimation, this simplification implies that r provides no information to help estimate x &++1 definition of shared and private states.
for all k; this understanding is consistent with
3.2 Prioritized generalized-linen dynamical modeling (PGLDM)
Our method, which we term Prioritized Generalized- Dynamical Modeling (PGLDM), uses a multi-staged learning approach to model primary generalized-lines time-series while prioritizing identification of the dynamics shared with a secondary time-series. Note, "primary" refers to the data source whose modeling is of primary interest and that can optionally be used to predict the

---

secondary data source. For example, within the context of decoding continuous behaviors from discrete population spiking activity. Poisson observations (i.e., neural activity) are the primary time- series whereas Gaussian observations (e.g., kinematics) are the secondary time-series. During stage 1, shared dynamics are learned using both observations. In stage 2, any private dynamics in the primary time-series are optionally learned. This two-staged approach allows prioritized learning of shared dynamics in the sense that latent states will be dedicated to explaining non-shared dynamics in the primary time-series only if there are enough latent states to explain the shared dynamics. Finally, an optional stage 3 allows identification of the dynamics private to the secondary time-series.
Below we will outline the first two stages of our algorithm and leave the optional stage 3 to appendices A crucial component of our method is a new covariance- based SSID algorithm for Identifying shared and private dynamics between two observation time-series using their first and second moments only. This covariance-baser algorithm what enables our method be applicable to generalized- observations, which we will expand further in section 3.2.3 We first present this new covariance-basee SSID algorithm, or equivalently PGLDM for linear state-space models with continuous Gaussian observations (the first three lines of equation (6)). before showing support for generalized-linee time-series broadly. The derivation of PGLDM provided in appendix A.I
3.2.1 Stage 1: shared dynamics
In the first stage, our algorithm identifies the parameter set corresponding the shared dynamical subspace, (A11, C( C! b, d), given hyperparameter nh and using both 2₃ and rg, both of which are observable in the Gaussian case (i.e., equation (1 or equivalently the first three lines of equation (6)) To do this, we first construct a Hankel matrix between future observations of the secondary process and past observations of the primary process
Azri
Agri-1 Aar:
Agn
Agr--1
Agra
H Cov(
(8)
AZra-11
Azrai-2
Agr.
(). Although equation uses the same horizon for both observations, in practice we implement the method for a more general version with distinet horizon values i, for the primary observations and ig for the secondary observations, resulting in Her E Risen, *i,* This allows users to independently specify the horizons for the two observations, which can improve modeling accuracy especially if the two observations have very different dimensionalities (see section 4.2 and appendix A. After constructing H we decompose using SVD and keep the top n1 singular values and their corresponding singular vectors
with 2, := [z
-- and rр defined as in equation
C₁ [A
CzAn
H r (2
AnG GOD
LC,
where nı is the user-specified dimensionality of the shared latent states x (1) I2 denotes the observ-
ability matrix for the secondary observations, and A (1) denotes the controllability matrix associated
with the shared latent states (defined in equations 3 and (16) in the appendix). At this point, we
extract C by reading off the first n, rows of Tg. To extract C we first form Hp per equation (2)
and extract the observability matrix for associated with the shared latent dynamics, ((1) by right multiplying H. with the pseudoinverse of A (1)
C c
An
HAA (Uf = r
c
A
We then read (( from the first n, lines of T (defined as in equation (0) the appendix). The baseline parameters b and d are empirically computed as the means of F3 and 2g, respectively.
5

---

Lastly, to learn the shared dynamics summarized by the parameter A11, we solve the optimization
problem A (1) A where (1) and A denote A from which n, columns have been removed from the right left, respectively. The closed-form least-squares solution for this problem is
An = A (1) )t This concludes the learning of the desired parameters (A11, C c( b, d),
given hyperparameter 11
3.2.2 Stage 2: private dynamics in primary process
After learning the shared dynamics, our algorithm can learn the dynamics private to the primary process that were not captured by x²" Specifically, learn the following parameters from equation
([A21 A22] c(2)). ) with hyperparameter n2 determining the unshared latent dimensionality of To do so, we first compute a "residual" Hankel matrix, H1 using rf") and A (1) from stage and decompose it using SVD, keeping the first 122 singular values and vectors
H = H, - r (1) r (2)
(10)
With 2) which corresponds to the first n, rows of ? we construct C, = c
(2) We
then A to form the controllability matrix A (1.2) the concatenation of A (1) and A (2)
(derivation in appendix A.I:
[A (1)
A (1.2) (A.22) ((122)
A G(12)
G(I = A (2)
where A (1,2) refers to the upper left block of the dynamics matrix A that corresponds to the
latent states x(1) and (2) Given A (1,2), we extract (A21 Azz] by solving the problem A (2)
(1,2)
[A21
A22
where
(2)
Az2 A022 G(1
[A21
Aaal GO:2)),
= (1.2) G(22) (
G
An
0
Concatenating the sub-blocks together, A (1.2) =
Thus, given hyperparameters nı and 72, we now have all model parameters associated with the shared dynamics and dynamics private to the primary signal: (A()) Cr, ((1) b, d). The remaining model parameters (A13, C2 are
A21
A22
learned in stage (appendices 1.4 and A.9
3.2.3 Supporting generalized-linear processes
The covariance-biic based SSID algorithm that we have just derived is what enables our framework to be broadly applicable to generalized-linen time-series data. As discussed in section 2.2 the variables that are linearly related to latent states x are unobservable in generalized- models.
However, because the algorithm outlined in sections 3.2.1 3.2.2 only relies empirical covariances and cross-covariances of the two observation time-series, we can support generalized- processes by using moment-conversiins (e.g., section 2.2 21 24 When computationally tractable moment conversion equations exist, we can compute both a Hankel matrix H,. as described in section 2.2 and also a cross-term Hankel matrix Her- For example, for the scenario wherein the Poisson observations constitute the primary process and Gaussian observations the secondary process (i.e., the first four lines of equation (6), we can compute a moment conversion to estimate joint moments of zı and r& from the joint moments of the observed signals 2₁ and yı with the following equation (derived using the conditional statistical properties, see appendix A.1.5)
(11) where, similar to equation 5). m and n correspond to indices of the first and second moments of the observation vectors z and rp (or yp), respectively. As another example, both generalized- observations (e.g., and t. from equation (6)) are Poisson distributed, the joint moments can be computed as
A Sm Tyn
Cov( Vp..) Pyre
Az
In(Cov(t/n.Vr. Mt,,.. Myp..) In My...)
(12)
For both of these scenarios the baseline log-rates are learned as in PLDSIID (section 2.2 Thus, our novel multi-staged covariance-bart based SSID learning algorithm enables identification shared vs. private dynamics across various generalized-linee processes.

---

Table 1: Shared mode identification error (log10, 10, i.., -2 means 1%) at the shared latent dimensionality
(ng = n1). indicates that a method (row) does not support the primary observation model (column).
Primary time-series (Tx or yı) / Secondary time-series (Zg or t²)
Method Name
Gaussian/Gaus.
Poisson/Gaus.
Pois./Pois.
BernoulliGGas.
PGLDM (stage 1)
-2.757 0.070
-2.707 + 0.091 -1.969 + 0.079
-2.664 + 0.072
Laplace-ENt [26]
-1.320 + 0.091
-1.083 + 0.119 X
-1.08 + 0.110 X
-1.027 + 0.067 -1,203 + 0.112
PSID (stage 1) [14
-2.985 + 0.102
x
Covariance SSID 19
-1,467 + 0.080
X
X
X
PLDSID 21
X
-1,319 + 0.132
X
bestLDS 24]
X
X
-1.209 + 0.117
4 Experimental Results
4.1 Shared dynamics are accurately identified in generalized- linear simulations
To evaluate how well our method identified the shared dynamics between two generalized-lines time-series, we simulated observations from random dynamical models as per equation (6). All state and observation dimensions were randomly selected, and the corresponding system parameters were randomly generated to simulate stable and slow-decay deccyeinnnn dynamics (see appendix A.7.I for details). In our first experiment, we evaluated how well the shared dynamical subspace could be identified when learning models at the true shared dimensionality. We performed this analysis for four combinations of generalized-linen observation pairs (r': or Yı and 23 or tx respectively):: (1) Gaussian/Gaussiaan (2) Poisson/Gaussion (3) Poisson/Poisson, and (4) BernouliiGaussian, Within each configuration, we compared models learned by our method against models learned with either Laplace-EM
maaiiiiiiimmmiiiiiiiiiinnn [26] or a SSID algorithm with the appropriate observation case we also compare against PSID [14,, an SSID algorithm that preferentially learns the shared dynamics between two Gaussian time-series. PGLDM (our method) and PSID were trained using both the primary and secondary time-series, whereas all other methods used only the primary time-series as they only model a single data source. We evaluated identification of shared dynamics by computing the normalized eigenvalue error between ground truth shared modes (i.c., eigenvalues of A11 in equation (7) and the identified modes (i.c., the learned A11 for PGIDMM/SIDD or A for the other baselines); the results of this analysis for 20 systems per configuration in Table For almost all conditions PGLDM more accurately identified the shared dynamics.
distribution. For the Gaussian/Gausse
appendix A for evaluation details. We report
In our second simulation experiment, we studied the effect of latent state dimension on learning. We generated 16 systems with fixed dimensions for shared and private latent states given by ni 4, n2 = 12, and rig- = 4, accordingly. We swept the learned latent state dimension from 1 to the true dimensionality of the primary observation time-series n1 + n2 16, with the dimensionality of shared dynamics set to minncuurent n2, n1). .. We found that our method accurately identified the shared modes with the minimal latent state dimension of 4; in contrast, PLDSID and Laplace-EN EM did not reach such high accuracy even when using higher latent state dimensions (figure (p) In these simulations also evaluated the predictive power (i.e., correlation coefficient, CC) of the model when using discrete Poisson observations to predict continuous Gaussian observations in a held-out test set (see appendix (A.8.2). This second metric allowed us to test our hypothesis that PGLDM's explicit modeling of the shared subspace improved decoding of Gaussian observations from Poisson observations compared with our baselines PLDSID (21) and Laplace-EM (26). We observed that our method achieved higher decoding performance in low-dimensional regimes, even when using as few as 4 latent states, whereas PLDSID required much larger latent state dimensions (around 12) to reach comparable performance (figure [b). We also evaluated Poisson self -reeicediition using area under the receiver operating characteristic curve, AUC (figure [p)). With the inclusion of stage 2 and sufficient model capacity, models learned by PGLDM were able to achieve comparable performance in self- eelf-prediction as compared our baselines.
7

---

(a)
(b)
0.78
Gaussian Prediction
Correlation Coefficient
0.7
for Poisson Self- Self-prediction
Oracle model performance
AUC
Oracl model performance
0.5
0.70
number shared latent states
True number of shared latent states
0.3
0.62
10
12
14
16
2
10
12
14
16
Number of Latent States (nx)
Number of Latent States (nx)
C
PGLDM
(d)
Normalized Shared Mode
Identification Error (log10)
-0.5
PLDSID
Polsson Laplace-EM
-2.0
Shared modes
Y-only modes
Z-only modes
PGLDM
-3.5
X
bestLDS
Bernoulli Laplace-EM
4
10
16
Number Latent States (nx)
Figure 1: In simulations, PGLDM better predicts Gaussian observations from Poisson observations, especially in low-dimensional regimes. Solid traces show the mean and shaded areas denote the standard error of the mean, (s.e.m.) for each condition. (a-b) Predictive power as function of latent state dimensionality for all learned models compared against oracle model. i.c., model with the ground-trume parameters. Left panel (a) shows prediction CC for the Gaussian observations and right panel (b) Poisson self-pre self-prediction AUC. (c) The normalized identification error of the shared dynamical modes (in log10 scale) as function of latent dimensionality. (d) Mode identification with models of size
accurately learns the shared dynamical modes and
sample Bernouli/Gauussia system with true dimensions n2 n2 6, n3
4.2 Modeling shared dynamics improves motor decoding from population spiking activity
As demonstration on real data, used PGLDM model the shared dynamics between discrete Poisson population neural spiking activity and continuous Gaussian arm movements in a publicly available NHP dataset from the Sabes lab [17]. The dataset is of a NHP moving 2D-cursor in virtual reality euvironment based on fingertip position. We use the 2D cursor position and velocity as the continuous observations z. For all methods we used 50ms binned multi-unit spike counts for the discrete observations We evaluated decoding performance of learned models using five-fold cross validation across six recording sessions (see appendix A.7.2 valss-valiiation details). For PGLDM, we use the shared dynamics dimensionality of n1 min ccurrent n,,, 8), i.e., maximum n1 of 8, because behavior decoding using stage roughly plateaued this dimension.
Compared with PLDSID and Laplace-EM, our method learned models that led to better behavioral decoding all latent state dimensions, including the maximum latent state dimension (figure 2 This result suggests that our method better learns the shared dynamics between Poisson spiking and continuous movement observations due to its ability to dissociate shared vs. private latent states. Interestingly, despite the focus on learning the shared latent states in the first stage, PGLDM was also able to extract the private latent states in the Poisson observations because of its second stage. This stage led to improved neural self- predictiioo AUC, while still maintaining the more accurate behavioral decoding (figure 2b-c). Indeed, even with the inclusion of just two additional latent states to model private Poisson dynamics (nz = 2, nl 10), neural self prediction was approaching that of models learned by PLDSID (figure (D). Finally, given its analytical nature, PGLDM required substantially lower training time compared with Laplace-EM EM (see appendix Table 2).

---

Kinematics prediction CC
Neural self-prediction AUC
0.67
0.65
PGLDM
0.25
PLDSID
0.63
Poisson Laplace-EM
0.15
n1=8
0.61
0.59
8
12
10
12
Number of Latent States (nx)
Number of Latent States (nx)
(c)
Kinematics prediction
0.70
Neural self seff-prediction
0.70
Cursor
Cursor Y-position
x- x-position
CC
AUC
x-velocity M Apbotpolt
Cursor
Cursor y-velocity
nx=12
nx=12
10 seconds
Figure In NHP data, PGLDM improves movement decoding from Poisson population spiking activity. (a) Solid traces show the average cross-validated kinematic prediction CC (shaded areas denote the s.c.m.) for models of different latent dimensions learned by PGLDM, PLDSID, and Laplace-EM. (b) Same as (a) but visualizing one-step ahead neural self- self-prediction AUC. (c) Kinematic prediction CC and neural self- self-prediction AUC for models latent dimensionality nt 12. Asterisks indicate statistical significance (Wilcoxon signed-rank test) with *; 0.05 and ***, p 0.0005. (d) Example decoding of cursor (x.y) position and velocity from test data.
4.3 PGLDM models better decode spiking activity of one visual area from another
As second demonstration on real data but with different combination of observation distributions (Poisson/Poisson), we used PGLDM to decode neural population spiking activity in one visual area from another. a publicly available dataset from Zandvakili and Kohn 16. 18 simultaneous V1/V2 population recordings were performed in anaesthetized NHPs as they were presented visual stimuli. We used five-fold cross validation to evaluate learned model performance decoding V1 activity from V2 activity and in V2 self- oeff-rrediction We again compare with PLDSID and Laplace- EM. For all learning algorithms we tested four latent state dimensions such that nu = ni + 722 6, 8}. For PGLDM we used the first two stages, setting the shared dynamics dimensionality to mj = min (current nz 4). Similar to results figure 2 we chose maximum 72 of 4 because decoding roughly plateaued at this dimension. The conclusions were consistent with those in figure 2 modeling the shared vs. private dynamics by PGLDM allowed for better decoding of V1 activity while maintaining comparable self-prediction of V2 activity. Analysis details in appendix A.7.3 We present the results of the complementary analysis (i.e., predicting V1 from V2) in appendix A.12
4.4 Limitations
PGLDM, similar other SSID methods, uses time-invariant model which may not be suitable if the data exhibits non-stationarity, such as chronic neural recordings. To handle non-stationarities. would need to either intermittently refit the model after a predetermined duration of time, or develop adaptive extensions 27 28). As an example of the latter, one can gradually update the model parameters by incorporating learning rate that weighs recent observations more heavily in the moment computations while gradually forgetting past observations [27 Moreover, as with other covariance- based SSID methods, PGLDM may be sensitive the accuracy of the empirical estimates of the first- and second-order moments. However, with increasing number of samples these empirical estimates will approach true statistical values, thereby improving overall performance, as seen in appendix figure 4 Due to errors in the empirical estimates of the covariances, SSID methods may also occasionally learn unstable dynamics (see appendix A.5). Future work may address this by

---

(6)
decoding AUC
0.60
V2 self-prediction AUC
0.550
ni=4
0.58
PGLDM
0.540
PLOSID
0.90
Laplace-EM
Polsson
0.530
0.54
n1
0.520
Number Latent States (nx)
(nx)
(c)
(d)
0.80
V2 022sfff-pee solf-prediction AUC
0.65
0.60
0.55
0.50
0.50-
0.40
0.45
nx=8
Figure 3: In NHP data, PGLDM improves V1 decoding from V2 population spiking activity while maintaining comparable V2 self-prediction performance. (a) Average cross cross-validees validated V1 decoding AUC (shaded areas denote the s.c.m.) for models different latent dimensions. (b) Same as (a) but visualizing V2 one-step ahead self-prediction AUC. (c) VI decoding AUC at n, = 8. Whiskers correspond to s.e.m. Scatter points are individual trials. (d) Same as (c) but for V2 self-prediction.
incorporating techniques from control theory, such as mode stabilization and covariance matching [29] Finally, although the GLDMs that PGLDM learns widely used (e.g.. in neuroscience), such models may not be suitable for time-series with nonlinearly evolving states. We did not compare our method against nonlinear deep learning methods, such recurrent neural networks and transformers 38] because the goals of these two modeling approaches are different. While nonlinear deep learning methods are typically used boost overall decoding performance, GLDMS used for their interpretability and utility in scientific investigations and in real-time, compatationally-effiii engineering applications (e.g., brain-computer interfaces).
32
5
Discussion
We developed PGLDM, a novel analytical multi-staged covariance-basse SSID algorithm for mod- eling two generalized- processes while also dissociating shared from private dynamics. In simulations we demonstrate that our method successfully achieves this capability agnostic to the generalized- observation distribution. As a result, our approach more accurately models system dynamics compared to several commonly-used GLDM variants and their corresponding learning algorithms. We also demonstrate our method's 's applicability to real data by modeling two distinct NHP datasets recorded under different contexts and from different brain regions. In both simulations and real data, PGLDM's ability to dissociate shared from private dynamics improved decoding of sec- ondary time-series from primary time-series despite using lower-dine divee---iensionall latent states. Further, although here we specifically focused on modeling Gaussian, Poisson, and Bernoulli observations, algorithm can be extended to alternate distributions described with generalized -anea models or to other link functions than the ones used here, as long as there exists a corresponding computationally tractable moment conversion equation. This is possible, if closed-fone equation exists, because the covariance- based approach of PGLDM only requires the second-order moments after moment conversion (equations (2). (5). 88 (II. (I2): ): such, in these scenarios the moment conversion algorithm can be modified for the desired link function and/or generalized-linear observation model 21 24
Beyond neuroscience, due the high-prevalence of GLDMs across various application domains, our method may be a useful tool for modeling the shared and private dynamics of joint generalized-linee processes with distinct observation distributions.
10

---

Acknowledgments and Disclosure of Funding
This work was supported by Office of Naval Research (ONR) YIP grant N00014-19-1-2128, Army Research Office (ARO) MURI grant 1911NF-16-1-0366 and National Institutes of Health (NIH) grants DP2MH126378, RO1MH123770, and Brain Initiative R61 R61MH135407. M.M.S. is an inventor on University of Southern California's 's patents or patent applications related to decoding and closed- loop control approaches, and is a consultant for Paradromics Inc. Authors thank Han-Lin Hsieh, Christian Song, and Parima Ahmadipouranari for helpful discussions.
References
[1] Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu. A look state-space models for neural data. Journal of computational neuroscience, 29(1-2):107-166, August 2010. ISSN 1573-6873 0929-5313. doi: 10. 10.000//00227-00--00117
[2] J. H. Macke, L. Buesing, and M. Sahani. Estimating state and parameters in state space models of spike trains, page 137-159. Cambridge University Press, 2015. doi: 10.1017/CB0978 1339941433.077.
[3] Scott Linderman, Ryan P Adams, and Jonathan W Pillow. Bayesian latent structure discovery from multi-ncuron recordings. In Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
[4] Wilson Truccolo, Uri T. Eden, Matthew R. Fellows, John P. Donoghue, and Emery N. Brown. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of Neurophysiology, 93(2):1074-1089 2005. doi: 10. (0.1155//n.00667.220 20047.2044 PMID: 15356183.
[5] Jonathan W. Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M. Litke, J. Chichilnisky, and Eero P. Simoncelli. Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature, 454(7207))9955 999,,999,.44 August 2008. ISSN 1476-4687. doi: 10.1008/natur00711 140.
[6] Vernon Lawhern, Wei Wu, Nicholas Hatsopoulos, and Liam Paninski. Population decoding of motor cortical activity using a generalized linear model with hidden states. Jowrnl of Neuroscience Methods, (89(2)2677-280, June 2010. ISSN 0165-0270. doi: 10. jn.euuhhjjee www.jnemehh2111.0...000. 2010.03.024.
[7] Hung- Yun Lu, Elizabeth S Lorene, Hanlin Zhu, Justin Kilmarx. James Sulzer, Chong Xie, Philippe N Tobler, Andrew J Watrous, Amy L Orsborn, Jarrod Lewis-Peacods and Samantha R Santacruz. Multi- scale neural decoding and analysis. Journal of Neural Engineering, 18(4):045013, Aug 2021. doi: 10. 1088/17411-2522411
[8] Yuxiao Yang. Shaoyu Qiao, Omid G. Sani, J. Isaac Sedillo, Breonna Ferrentino, Bijan Pesaran, and Maryam M. Shanechi. Modelling and prediction the dynamic respouses of large-scale brain networks during direct electrical stimulation. Nature Biomedical Engineering, 5(4):324-344 345, April 2021. ISSN 2157-846X. doi: 10. 0000//41551-00000000000
[9] Edoardo Balzani, Jcan-Paul Nocl, Pedro Herrero-Vide Vidal, Dora Angelaki, and Cristina Savin. probabilistic framework for task-aligned intra- and anter-area neural manifold estimation. In The Eleventh International Conference on Learning Representations, 2023.
[10] Jonathan C. Kao, Paul Nuyujukian, Stephen I. Ryu, Mark M. Churchland, John P. Cunningham, and Krishna V. Shenoy, Single-trial dynamics of motor cortex and their applications brain-machine interfaces. Nature Communications, 6:7759, July 2015. ISSN 2041-1723 doi: 18.08/m8mmmm.8.5..
[11] Mehdi Aghagolzadeh and Wilson Truccolo. Inference and decoding of motor cortex low-dimensional dynamics via latent state-space models. IEEE Transactions Neural Systems and Rehabilitation Engi- neering, 24(2):77--222, 2016. doi: 10. 1099/TNSRF..0115.2775522
[12] William Allen, Michael Z. Chen, Nandini Pichamoorthy, Rebecca H. Tien, Marius Pachitariu, Liqun Luo, and Karl Deisseroth. Thirst regulates motivated behavior through modulation brainwide neural population dynamics. Science, 364(6437)
April 2019. doi: 10.1 1 /26/science.aa 3932.
[13] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Mattco Carandini, and Ken- neth D. Harris. Spontaneous behaviors drive multidimensional. brainwide activity. Science, 364(6377); 789v7 April 2019. ISSN 0036-8075, 1095-9203. doi: 10. 12./

---

[14] Omid G. Sani, Hamidreza Abbaspourazad, Yan T. Wong, Bijan Pesaran, and Maryam M. Shanechi. Modeling behaviorally relevant neural dynamics enabled preferential subspace identification. Nature Neuroscience,
January 2021. ISSN 1746-1726. doi: 10. 10.1038/s411990 003000000000.
[15] Evren Gokcen, Anna L. Jasper, João D. Semedo, Amia Zandvakili, Adam Kohn, Christian K. Machens, and Byron Yu. Disentangling the flow signals between populations of neurons. Nature Computational Science, (18)):512--55 August 2022. ISSN 2662-8457. doi: 10. (0888/435880022-002222
[16] Amin Zandvakili and Adam Kohn. Coordinated neuronal activity enhances corticocortical communication. Newron, 87())822 839, August 2015, doi:
026..
[17] Joseph E. O'Doherty. Mariana M. B. Cardoso, Joseph G. Makin, and Philip N. Sabes. Nonhuman Primate Reaching with Multichannel Sensorimotor Cortex Electrophysiology. May 2017. URL https
//doi.orgg 10. /281/zenodo 553333
[18] Amin Zandvakili and Adam Kohn. Simultaneous V1-V2 neuronal population recordings in anesthetized macaque monkeys, 2019. URL http:s/s///err.....rraaaaaaaaaas wttp://aaaaa.ss////////eeeeaeeee
[19] Peter Van Overschee and Bart De Moor. Subspuce Identificattinn for Linear Systems. Springer US, Boston, MA, 1996, ISBN 978-1-4613-8061-0 -0. doi: 10.1007/978-1-4400 1---1----11156555 0465-4.
[20] Tohru Katayama. Subspace Methods for System Identification. Springer London, 2005. doi: 10.1007/
1- 1-84628-155 158-x.
[21] Lars Buesing, Jakob H Macke, and Mancesh Sahani. Spectral learning of linear dynamics from generalised- linear observations with application to neural population data. Advances in Neural Information Processing Systems, 25:9, 2012.
[22] Parima Ahmadipour, Omid G. Sani, Bijan Pesaran, and Maryam M. Shanechi. Multimodal subspace identification for modeling discrcte-continunt spiking and field potential population activity. Journal of Neural Engineering, 2023.
[23] Parsa Vahidi, Omid Sani. and Maryam M Shanechi. Modeling and dissociation of intrinsic and input- driven neural population dynamics underlying behavior. Proceedings of the National Academy of Sciences. 2024.
121(7)xe221288712)
[24] Iris Stone, Yotam Sagiv, Memming Park, and Jonathan W. Pillow. Spectral learning of bernoulli linear dynamical systems models for docision-making Transactions on Machine Learning Research, 2023, ISSN 2835-8856.
[25] Anne C. Smith and Emery N. Brown, Estimating a state-space model from point process observations. Newral Computation, 15(5):965 991, May 2003. ISSN 0899-7667. doi: 10.1 (00166228899666666622222222
[26] Scott Linderman, Benjamin Antin, David Zoltowski, and Joshua Glaser. SSM: Bayesian Learning and Inference for State Space Models, October 2020. UtLhhttss:///// //:::://hthttt /htggttttuuuuae. com/ hemm//inddrmanlllll sam
[27] Parima Ahmadipour, Yuxiao Yang. Edward F. Chang, and Maryam M. Shanechi. Adaptive tracking of human ECoG network dynamics. Journal of Neural Engineering, 18:016011, 2 2021. ISSN 1741-2552 doi:
[28] Yuxiao Yang, Parima Ahmadipour, and Maryam M. Shancchi. Adaptive latent state modeling of brain network dynamics with real- time learning rate optimization. Journal of Neural Engineering, 18(3)0030133, March 2021. ISSN 1741-2552, doi: 10.108///741-2255/aeeed
[29] J. M. Maciejowski. Guaranteed stability with subspace methods. Systems & Control Letters, 26(2): 153-156, 1995. ISSN 0167-6911. doi: https:://oo..or//1 0ttp://doo..o///00000010000661 1016/0167-69 (19(111090000))))
[30] Anders Lindquist and Giorgio Picci. Canonical correlation analysis. approximate covariance extension, and identification of stationary time series. Automatica, 32(5):709 -33, 1996. ISSN 0005-1098 doi:
https://doi.or//10.. 10t//111.11111111 1016/005-- 10699.
[31] Christopher Byrnes, Sergei V. Gusev, and Anders Lindquist. convex optimization approach to the rational covariance extension problem. SIAM Journal on Control and Optimization, 37(1):211 1-29,, 1998. doi: 10.1 1137/0066301999732
[32] Brien Alkire and Lieven Vandenberghe. Convex optimization problems involving finite autocorrelation sequences. Mathematical Progranmingg 93(3):331-359, December 2002. ISSN 1436-4646. doi: 10. 1001/ s10107-002-0334xx
12

---

[33] Chethan Pandarinath, Daniel J. O' Oheaa Jasmine Collins. Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman. Stephen Ryu, Leigh R. Hochberg, Jaimie M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential auto-cncoders. Nature Methods, 2018, doi: 10. 10.1038//11922 092-111-01911
[34] Cole Hurwitz, Akash Srivastava, Kai Xu, Justin Jude, Matthew Perich, Lee Miller, and Matthias Hennig. Targeted neural dynamical modeling. In Advances in Neural Information Processing Systems, volume 34, pages 29379-29900 293922. Curran Associates, Inc., 2021.
[35] Joel Ye and Chethan Pandarinath. Representation learning for neural population activity with Neural Data Transformers, Neurons, Behavior, Data analysis, and Theory, August 2021. ISSN 2690-2664, doi: 10.5 N.51162//001.223388
[36] Joel Yc, Jennifer Collinger, Leila Wchbe, and Robert Gaunt. Neural Data Transformer 2: Multi-context Pretraining for Neural Spiking Activity. In Advances in Neural Information Processing Systems, volume 36, pages 80352-80374. Curran Associates, Inc., 2023.
[37] Hamidreza Abbaspourazad, Fray Erturk, Bijan Pesaran, and Maryam M. Shanechi. Dynamical flexible inference of nonlinear latent factors and structures in neural population activity. Nature Biomedical Engineering, 800) January 2024, ISSN 2157-846X. doi: (0.1088/411511-003-01000 0016-11 1, Publisher: Nature Publishing Group.
[38] Omid G. Sani, Bijan Pesaran, and Maryam M. Shanechi. Dissociative and prioritized modeling of behaviorally relevant neural dynamics using recurrent neural networks. Nature Neuroscience, 27(10): 2033-2045, October 2024. ISSN 1546-1726. doi: 10. (0.1038644119-0024-0110 Publishing Group.
Publisher: Nature
[39] Steven Diamond and Stephen Boyd. CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1-5,, 2016.
[40] Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42-60, 2018.
[41] Mark M. Churchland, John P. Cunningham, Matthew T. Kaufman, Justin D. Foster, Paul Nuyujukian, Stephen I. Ryu, Krishna V. Shenoy, and Krishna V. Shenoy. Neural population dynamies during reaching. Nature, 487:51-56, 2012. ISSN 00280836. doi: 10. 10.1038/naturell 1129.
[42] Christian Y Song, Hin Hsieh, Bijan Pesaran, and Maryam M Shanechi. Modeling and inference methods for switching regime-dependent dynamical systems with multiscale neural observations. 2022. doi: 10. 111111881111-11 0088/174112555 255/4/49922////
[43] Daniel W. Moran and Andrew B. Schwartz. Motor cortical representation of speed and direction during reaching. Journal of Neurophysiology, 82(5))2666- 1999. doi: 10. 10.1152/jn.199 1999f/..99998988 0399.82.522677 PMID: 10561437.
[44] Shy Shoham, Liam M Paninski, Matthew R Fellows, Nicholas G Hatsopoulos, John P Donoghue, and Richard A Normann. Statistical encoding model for primary motor cortical brain machine interface. IEEE Transactions on Biomedical Engineering, 52, 2005. doi:
[45] Patrick N. Lawlor, Matthew G. Perich, Lee E. Miller, and Konrad P. Kording. Linear-nonlinear-time---- poisson models of neural activity. Journal of Computational Neuroscience, 45(3):173 191, December 2018. ISSN 1573-6873, doi: 10. 10.1007//1822-011 -01---66666-
[46] Shinsuke Koyama, Uri T. Eden, Emery N. Brown, and Robert E. Kass. Bayesian decoding of neural spike trains. Annals of the Institute of Statisticul Mathematics, 62(1):37-59, February 2010. ISSN 1572-9052. doi: 10. 10.1007/10466--00 084666300-02444 009-0229-x.
[47] Uri Eden, Loren M. Frank, Riccardo Barbieri, Victor Solo, and Emery N. Brown. Dynamic analysis of neural encoding by point process adaptive filtering. Neural Computation, 16(5):911-998 May 2004. ISSN 0899-7667. doi: 10. (22/08897664477113355559
[48] Skipper Seabold and Josef Perktold, statsmodels: Econometric and statistical modeling with python. In
9th Python in Science Conference, 2010.
[49] F. Pedregosa. G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer.
R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-lcam: Machine learning in python. Journal of Machine Learning Research, 12:2825-288 -230,, 201
13

---

[50] Joāo D. Semedo, Anna Jasper. Amin Zandvakili, Aravind Krishna, Amir Aschner, Christian K. Machens. Adam Kohm, and Byron M. Yu. Feedforward and feedback interactions between visual cortical areas differrntt population activity patterns. Nature Communications, 13(1):1099, March 2022. ISSN 2041-1723.
doi: 10. N0.103//41467-0222 022-7-222-2882--
14

---

A Appendix
A.1 Derivation
Here we provide the derivation for PGLDM (prioritized generalized- -inear dynamical modeling), covariance-based subspace identification algorithm that learns dynamical model of primary time-series while dissociating shared vs. private latents with secondary time-series. We define the following equivalent formulation for dynamical model (equation (6)) where the block structure delineates shared (x()) and private (x(2), x()) latent states
k+1
A 11
0
0
1++
(2)
A21
Ag2
0
(2)
We
0
Ass.
(3)
(3)
x x
(1)
Xx
rk
c
(1)
C
0]
X
+ V/:
(3)
%
(13)
(1)
X
(1)
0
(3)
x
+ E₁ + d
(3)
Yk Ir
Pyir g(rx))
tk
(tu h(zx))
with parameters and noise terms defined in sections 2.1 2.2 and 3.1
A.1.1 Standard Covariance- SSID
Before we present the derivation for PGLDM, we review a few steps in standard covariance-based SSID (section 2.1 that will help us in the derivation. First, it can be shown that the 7-th lag cross- covariance terms for r can be written in terms of model parameters as Ar. Cov (Fa4r: rx) = C,AT FEAT-1ca G, where G Cov(xx+1, rk). Using this relationship. the Hankel matrix, H,, can be expanded as 20
Ar.
Ar
A
H. Cov(ry,
A
A
(14)
C.A
C
C,G
C C
C,Ai-G
Second, using singular- value decomposition the above Hankel matrix H, be decomposed into observability, Tr, and controllability, A, matrices from which model parameters can be extracted
19 20
C,
C.A
SVD H, r
T,A=
AG G
(15)
C
We note that there exists canonical correlation analysis (CCA) version of this algorithm wherein H, left and right normalized by the square root matrices of the future and past observation covariance
15

---

matrices, respectively, prior the singular-value decomposition. Specifically I- H syp USVT, where バ Cov(rr, LL" and App Cov p. rp) MM The observability and controllability matrices can be recovered as H, = LÚ and A = 51/2/T respectively,
We refer readers to section 8.7 of Katayama (20] for more detail. We can similarly derive a CCA version of PGLDM, which we present next. This can be achieved by applying the appropriate normalizations to both H, and Hrr, the cross-term Hankel matrix defined below, prior to the singular-value decomposition. The resulting decomposition would then have be un-normalized to retrieve the appropriate observability and controllability matrices.
A.1.2 PGLDM: Stage derivation
In the first stage of our algorithm, our goal is to learn the model parameters that correspond to the shared dynamical subspace of z and r via the latent state (!); (A11. Cf C b, d). Before deriving the first stage, we first present few parameter definitions and simplifying assumptions. As
noted in section 3.1 G can partitionedd as G G ((2)" "
T
due to the block structure
of equation 13). Further, based on assumptions 3.2 and 3.3 we can set G 0. Leveraging the block structure of G' we can also define a simplified Block-partitions structure for the controllability
matrix A (equation 15 as
AG
G]
A11
0
0
An
0
0
G
GOT
A21
A22
Ar-GG
A21
A22
0
G
G(2
0
Aga
Ag3_
G
G3)
A (())
AnG
A (1)
AA21
A22
0) A
[A21
A22
G
G
(2)
(16)
33 A ((3)
G
(
(3)
A
A A (())
G
AA21
A22
0) A AG"
[A21
A2
(
2)
(
The block-partitin structure of A will play important role in the separation of learning between stages 1 and 2 and will enable the prioritized identification of shared dynamics A11-
To begin the derivation, it can be shown that the 7-th lag cross-covariance between and can be written in terms of model parameters as Agr, Cov (Zk+r, rR) = C.A -- GG. Using the block definition of G, we can simplify the cross-covariance term as
A 1
A Agr=
C (I
C
(3)
Az1
Az2
0
Aaa_
G
The Hankel matrix between future secondary observations and past primary observations can then be expanded as
Azr;
Azr,-1
Azri
H Cov(zf:Pp)
Agru-i
Agrac
Agr.
c (1) A A ((1)
C 2
ca (())
((l) 2 2 (
c A
2i-3 (())
c A (l-GG()
16

---

A singular-value value decomposition of Her yields the observability matrix for z (i.e., I2) and the controllability matrix A (1) associated with the shared dynamics
C
An
SVD H A1)
[Aii|GO)
An GD].
(17)
C2 ( A
At this point, C2 can be read off the first n, rows of Tzz. The shared latent dynamics matrix A1 can be learned by solving a least-squares problem based the controllability matrix A (1) (as introduced section 3.2.1)
(1)
Ani A
where
(18)
A = G
AnGoI
A
(AT
G
which has the following closed-form solution: A11
To extract c( we first note that the Hankel expansion in equation 15) can be rewritten with the block-structure of A in mind as (we omit partitions corresponding to latents xx() due to the decoupling)
L
T.
A
r. (UE
3
H, UEV*
VT)
(4)
(ww U
[
(19)
(U()))))
(I)) v
9(2)
(2)) v(2)7)
(1)
where equivalence (a) due to the dock-partition structure of A and equivalence (b) implicitly introduces a block structure to observability matrix Fr. where rfa and r correspond to the observability matrices associated with the shared and private latents, respectively. More formally,
C₁
C,A
(1)
c
r (1)
r
(20)
r
F(2)
C
where Т, denotes Tr from which the top rows have been removed. Taken together, H, can be viewed the sum of "shared" and "private" components (equation (19 Thus, we can compute r(")
as
H,A (it (T(2) A (1) (
(1) -!("))
(21)
where we have used the orthonormal property of right singular vectors V to conclude A (2) (1)
At this point, we can extract c(2) by reading the top πi, rows of ( Finally, b and d can both be learned directly from the data either by computing the empirical mean (when working with continuous Gaussian observations) during the moment transformation (section 2.2 This concludes the leam-
ing of all parameters associated with the shared dynamical subspace, i.e., (A C(1) ((() b. d).
A.1.3 PGLDM: Stage 2 derivation
In the second stage of our algorithm, our goal is to learn model parameters that describe the private dynamics of via the latent state x (2), ([A21 Aza] (2) ). To learn these parameters, first
extract the private component in equation 19 termed H. (2)2) by subtracting ( (1) from H,,
17

---

and decompose it via singular-value decomposition to get ( and A (we omit reference to private latents (3) due to the decoupling)
= H, - (1) (2) (2)
(22)
At this point, we take Cf2) as the top n,, rows of r and concatentate with C to complete C,, To complete the state dynamics matrix A, we refer back to the block-structure representation of the controllability matrix in equation (T6)
A (1) (2)
An
A11
G G
G G
A
A21
A22
A
Az
from which construct the following relationship
A (1)
A11
A (1)
(23)
A
A
A22
A (2)
where A and A are defined as in equation (8).
can further isolate the residual state transitions the solution to the following equation (taken from the second row of equation (3
(2) = A2
A22]
(2)
=
[A21
Agu
A
(24)
which has the following dlosed-form least-squares solution:
[A21
An2
A
A2)
The full state
A11
dynamics is the concatenation A
This concludes the learning of parameters for the
A
A2a
private dynamics in i.e., ([Ag1
A22].
(2)2)
A.1.4 PGLDM: Stage 3 derivation
Finally, in the third state we can learn the model parameters associated with the private dynamies in that is (Aga, ((3)), ). in an approach similar to stage 2. We first construct future-past Hankel matrix, Hz. associated with the secondary observation z. Using analysis similar equations (19) and (2T). we can show that H2 as the sum of "shared" and "private" components. the case of secondary observation z, we can explicitly show that its observability matrix, I assumes block form as
C₁ C.A
C C
33)
r₂
(a)
A
C
3
A33
r
O
r
(25)
C
(
A--1
C
1-1 33
where equivalence (a) is due to the block dookt-patiiionn definitions C2 and A in equation (13). forth, we omit the middle block of Tz, without any loss of generality, and simplify H, as
Hence-
[r ria
H2 UEVT (3
(UX) (EI/VTT)
Enyy
vaa"
(10))
U
0
(/)
0
23
[V()'
rp (U
r (U
(3))
(3)
(3)
v(3)")
(3) (3)
(26)
where equivalence (a) used the block structure of T, and equivalence (b) implicitly introduced block structure on the controllability matrix associated with observation z. In stage we had learned
18

---

T thus we can use an approach similar to equations (1
22 to extract the private component of
Hg. We first compute A the controllabilty matrix associated with the shared dynamics, as
T H2 - ( - + + (3) =
(27)
Then we can extract and decompose the private component of the Hankel matrix
H = H r - A sVD r A
(28)
(33) taken as the top ni rows of T thus completing Cz. To learn A33 we the same approach standard covariance-bate based SSID (section 2.1 and solve the problem F = p As3. where Fg and F3 denote Fz from which the top or bottom n2 rows have been removed, respectively. This optimization problem has the closed-form least-squares solution A33 T.Tz z. This concludes the learning of parameters for the private dynamics in z, i.e.,
As in the previous stages, C
(A33, Cx
A.1.5 Transformation of joint Gaussian and Poisson moments
In sections A.I.1. we demonstrated how all model parameters be extracted in three stages with prioritization, starting from the moments can be computed from generalized-lines Poisson/Gaussiar observations as and z, respectively. The joint second-moment Az5.,. Fya, can be analytically recovered from the computable moments of y (the Poisson observations) and (the Gaussian observations), using equation
of the and z. As an example of how these consider the specific case of observations,
as
Ax Cov(
where m and n indicate index-wise notation. Here we provide sketch of the proof. Without loss of generality, assume and are stationary with mean of 0 (e.g., demeaned during preprocessing). We can compute the covariance of any two elements j and le of vectors 2 and yp respectively as
Cov(2 ype)
E [21,NpL
E Z
[pe]]
25 PpL
E E Fpa]] = [ [Pp)] exp(xp=)]
where (a) is because 2 and You are independent when conditioned on latent log-rate ry. Next, we use the fact that zy and rp are jointly Gaussian random processes and, result, the mean of the conditional distribution, Elz,, [rp... equal Axy, tMx A rps (i.e., the linear least-square estimate of z using p.,)) The last step is compute the expectation
E [ A ppa exp(rpu) Azi,Foy Frr Myr. Cov(zzy. Vpu)
which, after rearranging terms, yields
Az Fun
Cov(Ezn: Vp..) Myy..
We note that the final equation is equivalent to a derivation provided by Buesing et al. 21 as their supplementary equation (6) compute cross-covarianees between Poisson observations and Gaussian imputs, instead of between joint Poisson and Gaussian observations (as was in our case). The remaining unimodal (i.e., Poisson-only) moment conversions that are required to compute H, are performed per equation (5) in section 2.2 For Aamouuli-Gaussian moment conversion equations, we refer the reader section 3.3 equations (5) and (6) in Stone al. (44)
A.1.6 Generalized cross-term Hankel matrix with different horizons per observation
For ease exposition, the derivation section A.1.2 was provided for a cross-term Hankel matrix Hyr that was formed with equal horizons for z and r as
Agr.
Agr Agr;
Agri
Azr
Agr
Zi
H = Cov(z),
Z
rp
Agrai
Agr.
Zz
r
19

---

In general, the rank of Hankel matrices formed from ideal data covariances can be shown to be
the same as the state dimension associated with it [99 20), i.c., n, = rank Hr,) per equation (17)
and n1 + 12 rank Hr) per equation (5) However, during system identification these Hankel matrices are formed from non-ideal empirical sample covariances and, as a result, are typically full
rank. Nevertheless, we expect the singular values associated with real dynamics (e.g., the first n singular values in H2r) to be larger than subsequent singular values that are due to noise. Indeed, the goal of the SVD applied to Hankel matrices, e.g., in equations (15) (17). and 222. is to remove noisy singular values and only keep the largest singular values that are most likely due to real dynamics.
Given that the Hankel matrices formed during system identification are typically full rank, their rank is determined based on their dimensions, i.e., rank Her) = min(i X n,, 2 x n2)) and rank (H,) 1 x Ry- Thus, the horizon parameter i that is used to form the Hankel matrix plays an important role in its final dimensions, rank, and, consequently, on the maximum number of non-zero singular values that can be preserved after applying SVD. This, in turn, determines the maximum state dimension that can be learned for the resulting model. Thus, to provide more flexibility over the state dimensions that can be learned in each stage of PGLDM, we generalize the Hankel matrix Her to support different horizon values for each of the observations, ig and ip, such that
A Zri
Arri Azr.,
A zr: 611
+
Azr., -4,2
Ziz
Har
with By =
Tp
Zei.
Agra.-1
A Agra.-2
The observation horizon i, is also used when forming the Hankel matrix Hr, per equation (1). The additional flexibility gained from having different horizon values can be especially critical in scenarios wherein the dimensionalities of and are very different, such as in the case of our NHP analysis in section 4.2 where n, = 4 and n, = 15. We select the final horizons i, and é, via an inner rross-validation based which values achieve the best decoding accuracy in the training data.
A.2
Assumptions and generalizability of the block structure formulation
Here we explain the assumptions underlying the proposed block formulation of A, as defined in equation (7). and discuss the model's generalizability. First, we assume that the private latents of y (i.e., x x) and of t (i.e., (3) do not drive the shared latents (i.e., )) so not to leak private dynamics into shared dynamies. Thus, we take A12 = and A13 = Second, we assume that the private dynamical latents of the two signals (i.e., (2) and x)) do not drive each other, and so we take A23 = O and A32 = O. Finally, the asymmetrical decision to set A21 * O and A31 = 0 was dictated by the asymmetrical roles played by the primary (i.e.. y) and secondary (i.c., t) time-series in our formulation. We initially made the assumption that A31 = 0 to simplify the derivation of our method because it allowed the optional third stage (i.e., extraction of the dynamics private to t) to operate independently of the first stage (i.e.. extraction of the shared dynamics). However, this assumption was still in line with our notion of shared and private dynamics; Aaı # O, then the primary time-series (y) would be informative (i.e., predictive) of the private dynamics in the secondary time-series (t) through the influence of the shared dynamics (1)) on future (3)
Although A21 fundamentally has similar information leakage impact in the reverse direction, chose to keep the more general form of A21 # O. Unlike the As1 case, a non-zero A21 does not couple the mtute learning. This is because we designate a primary time-series in our formulation (taken as y here) that acts as a predictor for the secondary time-series (taken as t here). Given this designation, we exclusively calculate all Hankel matrices using the primary/predictor time-series as the past observations. Specifically, we only compute Axr, =
Fss) (equation (8b) and not A Cov(r 24). From a derivation perspective, this means that Ag1 only appears in stage and not in stage 1, thereby never coupling the learning of private and shared dynamics. However, we clarify that if the user requires a strictly symmetrical definition of A wherein all latent states are fully decoupled from each other (i.c., a block biack-diagonal structure), our algorithm can, with a minimal modification to stage 2, cover the case where A21 = 0. To learn block diagonal A, one
would only need compute the solution A (2) Az2 A (2) instead of the original formulation
A (2) [A
A2e
A (1.2)
(see section 3.2.2 and equation 24 ).
20

---

Finally, taking A21 # O also has particular significance in the case when there are no private latents associated with t, that is when n3 O. In this case, taking A2I non-zero allows the alock-structured state-space form defined in equations (6) and 7 to be assumed without any loss of generality. Specifically, when fi3 = O, the latent states our model describe the primary time-series y with a subset also explaining the secondary time-series t. Formally, we define the true dimensionality of the shared states (denoted by n1) based on the rank of the observability matrix for the pair (A, C.). It can be shown using linear systems theory that an invertible linear transformation of the latent states (i.e., a similarity transformation) always exists that can place the ni dimensional latent subspace that is observable via t as the first few dimensions of the latent space, thus giving the block-structured formulation of equation
This can be seen by applying Theorem 3.8 from Subspace Methods for System Idenrification (20] to the first two lines of equation (6). As a result, the blocked formulation of equation 7 is equivalent to the formulation from (6) and we can aim to learn our model in the form of equation 7 without any loss of generality. Thus, our choice keep A21 o 0 was also in
part motivated by a desire to maintain a more general state-space form for the case when fg = O.
A.3 Noise statisties
Standard SSID algorithms (e.g., section 2.1 learn linear state-space models of the following form
X Ax + Wi
(29)
rk = C,X: + Vk
where state noise, Wk. and observation noise, Vk. are typically additive Gaussian noise and may have a non-zero instantaneous cross-covariance S = Cov(W&, Vx). SSID in general does not assume any
restrictions on the noise statistics. However, the Poisson observation model (equations (4) and (6) has no additive Gaussian noise for тą and instead exhibits Poisson noise in Y& when conditioned on
Tg. This means that Væ = O in equation (6) and thus = 0 and S O. Imposing these constraints is important for accurate parameter identification for Poisson observations, but was not previously addressed by Buesing et al. (21) Thus, we require our algorithm to find a complete parameter set O' that is close to the learned (A, Cr,, Cz b) by PGLDM and imposes the noise statistic constraints R= and S = 0. To do this, inspired by Ahmadipour et al. [22], we form and solve the following convex optimization problem to satisfy the noise statistics requirements
minimize Ax
+ [ I
(30)
such that Ax . O, Q F O, R(Ax) 0
where Ax == Cov(X,, x²) denotes the latent state covariance. Further, we enforce the following covariance relationships, derived from equation
( [9], as constraints
{ ((Ax) Ax
AT
R(Ax) = Aro C,
(31)
S(Ax) = G
AA.C
This approach has multiple benefits. First, it finds noise statistics that are consistent with the assumptions of the model (e.g., R = 0). Second, it enforces the validity of learned parameters, i.e., parameters corresponding to a valid positive semidefinite covariance sequence (see section 4.44) It also enables state prediction (see appendix A.8.2). Combining the previously found parameters and the matrix Q that corresponds to the minimizing solution Ax of equation (30), we have the full parameter set O' = (A, Cr, Cz, b, Q). We used Python's CVXPY package to solve the semidefinite programming problem defined in equation (30) 33 40). For of our comparisons against PLDSID, we learned the noise statistics associated with the method's 's identified parameters using this approach, but kceping the rest of the algorithm the same.
A.4 GLDM parameter equivalencies
Here we briefly discuss the equivalence of learnable parameters by SSID for linear dynamical models of the form in equation (T). In section 2.1 we stated that standard covariance-based SSID learned the parameters θ = Cr G. Arg). but not any of the noise statistics (e.g., covariances Q and R). Indeed, there also exist SSID approaches that instead learn the following parameter set (A, Cr, Q, R, S) [19 20). Both sets of parameters are valid and knowledge of either completely
21

---

defines the model in equation T.. Moreover, the two sets of model parameters can be related to each other using equations (1 19 the same relationships we use as constraints in the convex optimization problem defined in appendix A.3 for learning the noise statistics of the PLDS model (T).
A.5 Possibility of learning unstable modes in small data regimes
Subspace identification methods generally only converge to the correct system parameters asymp- totically (see figure 4 in appendix A.IIId A. as the empirically estimated covariances also converge to their true values [9] For finite samples, however, there will always be some error in the learned parameters. Although such errors are generally benign, extreme scenarios can result in unstable state dynamics, i.e., the identified A has at least one eigenvalue with magnitude larger than 1. We omitted any learned unstable models in our analyses, which was reflected in the reduced number of samples in the standard error of the mean (s.e.m). However, we rarely encountered unstable models; for example, in the results presented in figure 4 there were unstable models for training set sizes typical of neuroscience datasets (i.e., le5 or le6 training samples).
A.6 Possibility of accumulation of error in multi-staged learning
There are two aspects of the learning to consider. First, there is signal-to-noiss (SNR) consideration that is fundamental to learning methods general. In stage 1, our method requires high SNR in the cross-correlations between future secondary time-series observations and past primary time-series observations. most of the signal present in the primary time-series is not attributable to the shared latent states (i.e., the residual Hankel in equation (TO) dominates), then most of H 's singular values will be small and stage may have greater estimation error. the reverse situation holds and the shared latent states explain most of the signal in the primary time-series, then the residual Hankel,
H will mostly have small singular values, possibly resulting in estimation errors during stage 2. Inspection of the singular values prior to model parameter extraction can, however, help guide method usage. For example, if the singular values of H2r are small, then this may indicate that the two time-series do not have shared dynamics and only stage 2 of the method is needed (i.e., standard GLDM). the singular values of the residual Hankel are small, then this may indicate that almost all dynamics are shared and only stage 1 is needed.
second consideration is a numerical one that could result in an accumulation of errors in downstream stages. If during stage 1 the modes are identified inaccurately because there is minimal shared dynamics between the two time-series and/or the training sample size is too small, then this estimation error could impact the computation of the residual matrix (2) thereby introducing error in stage 2. Examining the singular values of Har may, however, also help avoid this situation. Although the multi-stage learning can lead to error accumulation in some situations, it comes with the benefit that our method has the ability more accurately identify the shared dynamics (when they exist) in stage 1, compared to existing GLDM methods. Moreover, Figures [b. 2p, and suggest that the impact of error accumulation is not severe and can be situation- -euation-dependentt for example, self-prediction performance of PGLDM reaches that of PLDSID in Db, reaches very close to it in 2p, and exceeds it bb.
A.7
Experimental details
Below we provide details on how simulations were generated and about our real data analyses.
A.7.1
Simulations
simulated generalized-ling linear observations from random models as per equation (6). For the simulations used to generate the results in Table [] we fixed the
For our synthetic data in section 4.1
number of shared and private latent states as n2 2, n2 6, and nı 4. We randomly selected
the observation dimensions with uniform probability from the following ranges: 10 < nt < 15
and either S n2 S 10. when the secondary observation was Gaussian, or 10 < n. < 15, when Poisson. For the simulations used to generate figure 2 we fixed the latent dimensions as nı n2 = 12, and ng = 4 and randomly sampled observation dimensions with uniform probability from: 20 < < 30 and 5 < n. < 10. Across all simulations we used these dimensions to generate random model parameters θ (A, Cr, Cz. b, d, Q). We constrained the complex eigenvalues (i.e.,
22

---

modes) of the state transition matrix A to have magnitudes uniformly distributed between [0.93, 0.99] and phases uniformly distributed between [0.019, 0.314]. These restrictions correspond to stable, slow-decaying systems with time-constants within [0.138, 0.995] seconds and frequencies within [0.3,5] Hz that are representative of various real time-series data, such as neural dynamics
For all Gaussian simulations we set the time-series mean (e.g., b or d) to O. Additive Gaussian noise covariances (e.g., Q) were randomly generated to be positive definite matrices. For observation time- series, the observation matrix C' was generated to hit a target signal-to-noise, defined as the variance associated with latent states normalized by observation noise variance (e.g.. [C2AxCe] //(e)))))))) Target SNR values were fixed at 10 for the simulations used in Table or were randomly generated as 10 with a uniformly distributed between [0, 2] for the simulations in figure | All Poisson simulations were generated on a 10 ms timescale with a baseline log rate (e.g., b) randomly selected within 10.5, 15] Hz. Observation matrix C' was scaled to achieve a desired per-dimension 95th percentile modulation depth of 10, where modulation depth is defined as exp(Cxx;) All Bernoulli observations were also generated on a 10ms timescale with a baseline set O. In this scenario, the observation matrix C' was scaled such that the corresponding Gaussian latent defined by Cx2 would have a standard deviation approximately equal to 1.
To simulate state noise covariance Q that adhered to the block-structure assumption defined in assumption 3.2. we leveraged the fact that latent (3) is completely decoupled from latents x(1) and x(2) in our model definition (as per equations 7 and (13). and assumptions 3.2 and 3.3). We chose to simulate the dynamics private the secondary observation time-series as separate 4-dimensional latent dynamical model. Specifically, we chose to generate one dynamical model corresponding to the latents x(2) and x²)), and a separate dynamical model corresponding to x(8). This approach, however, only works for Gaussian secondary observations. For the Piisson-Poisson case we instead generated an unconstrained Q. Despite this deviation from assumption [.2] our method's 's ability to
identify the shared dynamical modes was not impacted, as shown by the results in Table []]
Finally, for all simulation experiments we use horizon values of 10 for both observations by default. For the simulations used in Table
[ we generated 25600 training samples, and for the simulations used in figure I we generated 2e6 samples, splitting 50/50 into training and test data.
A.7.2 NHP Dataset 1: motor cortex recordings during reaching
All NHP analyses in figure 2 were performed on a public dataset released by the Sabes lab
17.. using the following sessions from monkey I: 201 60915/01, 20160916/01, 20160921/01, 20160927/04, 20160927/06, 20160930/02. We performed vross-validation using randomly-selected, non-overlapping subsets of 15 channels (n, = 15) binned at 50ms resolution within each session. We used a nested inner cross-validation to select hyperparameters per fold based on the prediction CC of kinematics in the training data. Hyperparameters in this context were discrete horizon ér. continuous horizon iz. and time lag, which specifies how much the neural time-series should be lagged to time-align with the corresponding behavioral time-series B3 43 44.. We swept i, values of 5 and 10 time bins, iz values of 10, 20, 22, 25, 28, and 30 time bins; and lag values of 0, 2, 5, 8, and 10 time bins. We removed channels that had average firing rates less than 0.5 Hz or greater than 100 Hz. Similar to Lawlor et al. [45), we also removed channels that were correlated with other channels using a correlation coefficient threshold of 0.4.
A.7.3 NHP Dataset 2: visual areas V1-V2 during stimulation
analyses in figure 3 were performed on a public dataset released by Zandvakili and Kohn using a randomly selected session 1071002p67. In the released dataset, each recording session consisted of repetitions of 1.28s of stimulus presentation (1 of 8 possible orientation gratings) followed by 1.58 of a blank screen. For our amalysis, we only considered neural activity during periods of stimulation binned at 20ms resolution. Further, we performed five-fold within each stimulation condition. For both observation time-series, we removed channels that had average firing rates less than 0.5 Hz or greater than 100 Hz. Because there were more VI units than V2 units, we performed cross- cross-validation using randomly-selected units equal to the number of available V2 units, that is n, = n2 = V22 units]. We present averaged results across folds, stimulation conditions, and oon-ovallapping subsets, which we collectively refer to as "trials" for a total of 120 trials. No hyperparameter optimization was performed for this analysis; discrete horizons for eSIDD-based methods were set to a fixed value of ir iz = 5.
16 18.
validvaadd modeling
non-overlapping subsets of V
23

---

A.8 Model evaluation
We evaluated learning using two different metrics: (1) shared dynamical mode identification for models with the true shared latent state dimensionality and (2) predictive power of inferred latent states.
A.8.1 Shared dynamical mode identification accuracy
To evaluate the learning of shared dynamics, we computed the normalized eigenvalue error as V| V/all2 / [ where Vune and 1a denote vectors containing the true (i.e., the eigenvalues of A1)) and learned shared eigenvalues, respectively. and || ||2 denotes the Euclidean 1.2-norm. To compute the first metric for our baselines that did not explicitly model shared dynamics, we needed to select the n modes identified from the primary time-series only that were the most representative of the secondary time-series. To do so, we first trained our baselines on the primary time-series observations and extracted the latent states. Then, we sorted these learned latent states based on their accuracy in predicting the secondary observations (appendix [A.822). We computed the eigenvalues associated with the top nı most predictive latent states, which considered the shared modes identified by our baselines. This was only necessary for configurations where the learned latent dimensionality was greater than the shared dimensionality (i.e.. nu n2 > n1). For configurations wherein learned n is smaller than true n1, we substituted missing modes with O prior to computing the normalized error.
A.8.2 Predictive power of inferred latent states
We also evaluate the predictive power of our learned models in two scenarios: (1) prediction of the secondary observation time-series from the primary time-series, common use-case in neuroscience 2 46). and (2) one-step ahead self- eeef-prediction of the primary time-series from itself. Using models learned from the training data, we constructed recursive Bayesian filters to estimate the latent states in a test dataset. Because our analyses in figures B used Poisson observations as the primary time-series, we chose use Poisson point-process filter (PPF) 47 for state estimation. (We note that using the PPF for state estimation is only possible if the learned noise statistics are valid, see appendix A..3) We denote the one-step ahead latent state prediction of x, using all samples of yı up to time k- ---1 - by X These state estimates can be used to predict both sets of (latent) observations as either CZ or C,
In order to compare the predictive power of our learned models with those of our baselines (PLDSID and Laplace-EM)) it was necessary to learn observation model (e.g., C₂2) for the secondary time- series. To do so, we first estimate the latent states in the training data using a PPF and then fit regression model (scikit-learn; statsmodels) from the latent states to the secondary observation For example, when the secondary observation, Zb. is a continuous Gaussian time-series, the parameter Cz was learned using ordinary least-squares such that C₂ = -X* (XX"),, where columns of Z and X contain 2 and & for all training timepoints k. To make all methods more comparable, we use the same approach to refit the secondary observation's 's model learned by PGLDM.
48 49].
For continuous Gaussian observations, we quantify the decoding performance using correlation coefficient (CC). For discrete Poisson observations we instead evaluate prediction using the area under the receiver operating characteristic curve (AUC). Since all Poisson predictions were made using the recursive Bayesian filter's ' estimates of the latent states, our goal was to validate if our model could accurately predict the occurrence point process events (versus no events) in a given time step when using all past observations of the primary time-series yı. For example, in the results shown in figure 2 we computed the probability of an event for the m-th dimension of y at time conditioned on all observations y1:k-1,
P > 01
Exs P(yT > 0 11 Y1=k-1)
(4)
Exe Yiuk- Ip(yT- > 0 xx)]
(b)
Exi yis-- [1 exp(exp(r")) xx]
(0) Exe LYsk-1 (exp(r?) xx]
(d)
mm [Hj + Aiim
exp
24

---

Table 2: Average learning runtime over 25 fits with 30484 training samples
Method Name
Runtime in seconds (# s.c.m)
PGLDM, 11 R2 = 8 (SSID/ / optimizzaiinn
0.269 + 0.008 0.080 = 0.005
PGLDM, n1 = 4, n2 = 8 (SSID/ / optimization)
0.253 + 0.005 0. 125 = 0.046
PLDSID, 12 8 (SSID/ / optimizaaion)
0.199 + 0.002 / 0.063 = 0.011
Laplace- Laplace-EM, f2 = 8 (100 iterations)
109.656 = 0.662
Laplace-ENa 72 = 8 (1 iteration)
1.097 = 0.007
where in (a) we simplify using y&'s conditional independence from the past Y1:k1++ in (b) we simplify based on yı Xx (exisso(exp((r)) in (c) use the Taylor series approximation of exp(exp(rx)) )) for small exp(rk). and (d) is simply the mean of a log-normal random variable. Note that F. C + b and A, C C where x is the current estimate for the state and Ax the estimate for the state- sateepredicction covariance. We can similarly compute AUC for the prediction of secondary Poisson observations, as in figure 3
A.9 Stage simulation results
We validated our algorithm's 's third stage for learning the dynamics private to the secondary observation z using simulations. We simulated 10 random models with the following dimensions n1 = n2 = aa = 4, n,, 20, n2 = 4, with the primary/secondary observation pair being Poisson/Gaussian. The Gaussian observations were set to have SNR of 10 (see appendix (pppendi(A...11// (A.7.I). We constrained the complex eigenvalues of the state transition matrix A to have magnitudes uniformly distributed between [0.93, 0.99] (as before) and phases uniformly distributed between [0.019, T.. Using a training set size of le6 samples, We then computed the normalized eignevalue error between the ground truth private modes (i.e., eigenvalues of Ass) and the identified private modes. We found an average normalized identification error of 1. 1.% for learning the dynamics associated with the latent states (3), For comparison, this is approximately equivalent to -2.05 on the y-axis of figure
used PGLDM with all three stages to learn all sets of shared and private modes.
A.10 Computation time details
We compared the computational runtime efficiency between our method PGLDM (using either stage 1 only or both stages 1 and 2), PLDSID, and Laplace-EM on one session of NHP data binned at 50ms resolution (section 4.2 and appendix (.7.2).. We repeatedly trained on 25 distinct time-series datasets and computed an average runtime for learning across all fits. Each dataset consisted of 6097-by-15 matrix (timestess-by-feature Gaussian observations. All methods learned latent models of dimension 8 (i.e., n2 = 8). When both stages of PGLDM were used, we fixed iterations (the default setting) and learning times are reported both as an average for 100 iterations as well as for single iteration. For PGLDM and PLDSID we separately report average running times for the SSID portion of the algorithm and for the convex optimization problem used to learn state noise, as outlined in appendix A.3 The results are presented in Table 2
of Poisson observations and 6097-by-4 matrix of
while keeping nt 8. Laplace-EM was run for 100
Most of the computational cost of our algorithm is involved in the matrix operations associated with 1) computing the necessary covariance/Hanker matrices, 2) performing the moment conversion, and 3) performing the SVD of the future-past Hankel matrices. To perform the moment conversion our method requires a covariance matrix for stacked future-past Poisson-Poisson observations (section 2.2 and a future-past Gaussian-Poisse of second-order covariances are computed using matrix multiplications which scale with the number of samples. As an example, we can consider the setup used for the computational cost analysis in Table wherein n, = 15 and i, 10 (horizon). The computed square Poisson-Poisson covariance matrix was of dimension 2 * n., 2, = 2 * 15 10 and was the result of a matrix multiplication between two matrices of dimension (2* 1110055-by-6000
Hankel matrix (section 3.2.1]. Both of these empirical estimates
where 6078- timesteps 2 i, Thus, this operation would scale linearly with the length of the training data. Similarly, the computational cost of this matrix multiplication scales linearly with feature dimension and horizon. The remaining operations (i.e., the SVD and the moment conversion itself) are functions of the latent-state dimension and the feature dimensions for each observation timeseries.
25

---

All running time analyses were performed on 2020 Macbook Pro using 2 GHz Quad-Core Intel Core i5 CPU with 16GB of 3733 MHz RAM.
A.11 Data Efficiency
We also investigated the data efficiency of our method by generating 20 random Poisson/Gaussian systems and studying the effect of training size on learning. The random systems used here were generated using the same procedure as the systems used in figure [ except for the latent dimensionalities of the primary signal. We instead uniformly sampled the total latent dimensionality of as < n1 + S 10 and the shared dimensionality as 1 S 11 S R2 + 122 we kept fixed size for the secondary observations private states R2 We used 1e3, 1e4, 1c5, 1e6 samples train models, and then tested them on lеб samples of independent held-out data (figure(4) We found that with increasing training set sizes the performance our method improved, both identification of shared dynamics and in overall predictive power. This is expected because empirically estimated covariances converge towards their true values with increasing training set sizes, thereby improving the overall performance of our algorithm.
(e)
0.80
(c)
Prediction Correlatio
Normalized Shared Mode
Iduntification Error (lea10)
Coefficiont (CC)
prodiction
0.78
0.76
-2.6
Oracie model performance
performance
+3.0
-3.4
0.74
3.0
4.0
5.0
6.0
3.0
4.0
5.0
6.0
5.0
4.0
Training sample size (log10)
Training sample size (log10)
Training sample size (log10)
Figure 4: PGLDM's overall performance improves with increasing number of samples, hitting peak performance with training set sizes of 1e5 samples.. Solid traces show the mean and the shaded areas denote the standard error of the mean (s.e.m.) for each condition. (a-b) Predictive power of learned models as function of training set size, with (a) depicting Gaussian observation decoding CC and (b) Poisson self-prediction AUC. (c) The normalized identification error of the shared dynamical modes (in log10 10 scale) as function of training size.
A.12 V1-V2 Results (time-series designation swap)
For completion we also present the complementary analysis the results in figure 3 wherein we take neural population spiking activity from V1 yı and activity from V2 as t.. The rest of the analysis kept the same as described in sections 4.3 and AA7. and the results are presented in figure 5 In this scenario, PGLDM does not make substantial difference in predicting V2 activity from VI activity when compared with PLDSID and Laplace-EM. Prior work using this same dataset has shown that the feedback direction, that is the V2 (past) to V1 (future) direction, generally exhibited higher correlations between the two time-series as compared to the feedforward direction [50]. We hottheesi that this finding might also explain the differences that in our results (i.c., the differences between figures and 5).
26

---

(b)
0.540
V1 self-prediction AUC
0.70
V2 decoding AUC
0.535
0.68
0.530
0.66
PGLDM
0.525
Poisson
0.62
n1=6
Laplace-EM
n1=6
0.60
0.515
Number of Latent States (nx)
Number of Latent States (nx)
(c)
(d)
0.60
V1 self-prediction AUC
0.85
V2 decoding AUC
0.58
95
0.54
0.52
0.50
0.48
0.46
nx=8
nx=8
Figure 5: All methods perform comparably in V2 decoding from V1 population spiking activity; PGLDM's V1 self-predictic performance improves with stage 2. (a) Average cross-validated V2 decoding AUC (shaded areas denote s.c.m.) for models of different latent dimensions. (b) Same (a) but for V1 one-step ahead self- celf-prediction AUC. (c) V2 decoding AUC at n2 8. Whiskers correspond to s.c.m. Scatter points are individual trials. (d) Same as (c) but for VI self-prediction.
27

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We state that aim in this work is to improve upon existing GLDMs, and their associated learning algorithms, that only consider a single generalized- time- series and therefore cannot dissociate shared private dynamics two generalized-linee time-series. We derive and present a novel algorithm with this capability and demonstrate in both simulations and real data, stated in the abstract.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We have a specific subsection called "Limitations" wherein we discuss the limitations of our model assumptions and methodology. We also have supplementary sec- tion that discusses the computational efficiency our method. Finally, our experimental details we discuss the number of simulated systems and amount (and type) of real data we have validated our method on.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best
28

---

judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We have full derivation of our method, with its associated assumptions, in our appendix. We provide a proof-sketch and outline of our method in the main manuscript.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We believe have provided all the necessary information regarding simula- tion configurations, open dataset locations, data preprocessing procedures, model hyperpa- rameters, and cross-validation settings needed reproduce the results. For the computational running time analysis we have also specified the hardware/compute resources used to run the analysis. Beyond reproducing analysis results, we have also outlined our analytical algorithm step-by-step in thorough detail to enable re-implementation by readers, if desired. Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
29

---

(c) If the contribution is a new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We have provided our current code with our submission and will provide public camera ready version upon acceptance. We also provide links to the open datasets used as well as to the codebase of one of our baselines.
Guidelines:
The answer NA means that paper does not include experiments requiring code. Please see the NeurlPS code and data submission guidelines (https://niis.cc/
cublic/guides//Cee htwpee/s/ooeeettbaiiiiiooooooooo
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //nips. wttpccccpppppu.....eeeee wwwicc/gguiss/eee /ooeeSe/e/oodbbmmiiiccea ssiinni asionPolicy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide details on experimental settings in the main text and in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
30

---

Answer: [Yes]
Justification: For all figure and table results we reported averages and their associated standard of measures. For figure 2 results also performed significance tests (Wilcoxon) and report statistically significant p-values.
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative еттог rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Most of our experimental results are agnostic to the compute resources used. For the computational efficiency results reported in the appendix do provide information on the compute resources used to run the experiments.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't it make into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http:////waaa.eecccpcpcccccccccc
Answer: [Yes]
Justification:
We have adhered to the conference Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
31

---

The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Although our primary intention with this work to develop novel SSID method from a foundations perspective, methods like ours have previously found utility in brain-computer interfaces and therefore work may have positive societal impacts by benefitting individuals with disability.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work does not make use of any models or data that could have a high risk of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
32

---

Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We provide citations and/or links to existing assets used in the manuscript including: scikit-learn, statsmodels, evxpy, the Python ssm package, bestLDS package on Github, and both publicly-available non-human primate datasets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: All model derivations and algorithm implementations described in detail in the manuscript. We have also submitted our code along with a README with instructions and tutorial that will be publicly released upon publication.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots. if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our work does not use crowdsourced experiments research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
33

---

According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/teview based on the requirements of your country institution) were obtained?
Answer: [NA]
Justification: Our work did not rely data collected from human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
34


Reviews:
[{"review_id": "lv9cJQCIKB", "reviewer": "Reviewer_Hfka", "summary": "This paper introduces a novel multi-step analytical subspace identification algorithm for Generalized-Linear Dynamical Models (GLDMs) to model shared and private dynamics within two time-series data sources. The proposed algorithm effectively decouples shared and private dynamics, demonstrating superior performance in simulations compared to existing methods. The algorithm's efficacy is validated using synthetic data and two non-human primate neural datasets, showcasing improved decoding accuracy of one time-series from the other.", "questions": "* I keep wondering what's the scientific meaning of modeling both the shared vs. private dynamics? Otherwise this would undermine the contribution of this whole work.\n* How does the computational complexity of the proposed multi-step algorithm compare with existing GLDM learning algorithms, especially in terms of runtime and scalability?", "limitations": "To my knowledge and understanding, there is no potential negative societal impact of this work. Other limitations please see “Weaknesses”.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* The paper is well-written and has a smooth and concrete flow.\n*  The paper presents a novel and effective approach to model shared and private dynamics with SSID theory.\n* The experimental section is sufficient.", "weaknesses": "* The multi-step nature of the algorithm may introduce significant computational overhead, potentially limiting its scalability.\n* The method assumes time-invariant dynamics, which may not hold for all neural data, potentially affecting its generalizability.\n* While the algorithm is validated on neural datasets, its applicability to other domains or types of data is not explored, which could limit its broader impact.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for finding our paper well-written, our approach novel and effective, and the experimental section sufficient. Below we address outstanding questions and comments inline.\n\n> The method assumes time-invariant dynamics, which may not hold for all neural data, potentially affecting its generalizability.\n\nWe thank the reviewer for raising this point and agree that the time-invariance assumption is a limitation. We included a discussion of this point in the Limitations section (line 277) and proposed two possible solutions: 1) intermittent model refitting, and 2) model adaptive extensions. For 1), one can refit the model after a predetermined duration of time, for example every day is typical in brain-computer interfaces. For 2), one can gradually update the model parameters by, for example, incorporating a learning rate parameter that weighs recent observations more heavily and gradually forgets past observations, similar to previous adaptive subspace methods [1]. We will revise the limitations section to include these discussions on possible extensions. We are also happy to include any other points the reviewer feels will be helpful.\n\n> I keep wondering what's the scientific meaning of modeling both the shared vs. private dynamics? Otherwise this would undermine the contribution of this whole work.\n\nWe thank the reviewer for their excellent question. We provide our explanation within the context of neuroscience but the same reasoning can apply to analogous situations from other application domains.\n\nModeling the shared dynamics is useful for studying how the brain encodes a particular behavioral/cognitive/affective process of interest, and for accurately decoding said process in brain-computer interfaces (BCIs). In general, the brain is simultaneously involved in several tasks or activities; for example, one may be moving their arm toward an object while also speaking and feeling excitement. As such, neural activity is a rich mixture model composed of multiple processes happening concurrently, thus necessitating the explicit dissociation of dynamics that are shared with a particular process of interest (e.g., movement). Another application, as reviewer m47S alluded to, is the development of BCIs. BCIs decode a behavioral/cognitive/affective process of interest from neural activity in real time. For example, BCIs can be used to restore lost motor function in paralysis, such as by decoding the patient’s movement intentions from neural activity to move a cursor on a computer screen. Dissociating and modeling the shared dynamics enables more accurate real-time decoding for BCIs as suggested by Fig. 2a, and is also more computationally efficient for real-time implementation due to the lower dimensionality required for the states, as suggested by Figs. 1a, 2a, and 3a.\n\nModeling the private dynamics separately, as an additional stage, is helpful for studying neural activity in its entirety (not just the shared component), for example to build a generative model of neural activity. Beyond generative modeling, this is also important for basic science investigations on the role of private neural dynamics during processes of interest. For example, some private dynamics that are not directly representing movement kinematics (e.g., shared with kinematics measures such as position and velocity) may be involved in time-keeping or in higher level cognitive planning during a movement task.\n\n> [From Weaknesses] The multi-step nature of the algorithm may introduce significant computational overhead, potentially limiting its scalability. [From Questions] How does the computational complexity of the proposed multi-step algorithm compare with existing GLDM learning algorithms, especially in terms of runtime and scalability?\n\nWe thank the reviewer for their great question. We have now expanded Table 2 to include computational time results for PLDSID and for PGLDM when using both stage 1 and stage 2 (vs. stage 1 only). To ensure consistency in hardware settings across all test conditions, we reran the entire experiment again and have updated the existing values in Table 2 to our latest results. We present the new table results here for convenience:\n\n|Method|Running time in seconds (mean &plusmn; STE)|\n|-|-|\n|PGLDM, $n_1=n_x=8$ (SSID / optimization)|0.269 &plusmn; 0.008 / 0.080 &plusmn; 0.005|\n|PGLDM, $n_1=4, n_x=8$ (SSID / optimization)|0.253 &plusmn; 0.005 / 0.125 &plusmn; 0.046|\n|PLDSID, $n_x=8$  (SSID / optimization)|**0.199 &plusmn; 0.002 / 0.063 &plusmn; 0.011**|\n|Laplace-EM, $n_x=8$ (100 iterations) |109.656 &plusmn; 0.662|\n|Laplace-EM, $n_x=8$ (1 iteration)|1.097 &plusmn; 0.007|\n\nOur method is more efficient than EM, which is iterative. Its run time is also only slightly higher than standard PLDSID (on the order of 10ths of a second). So in many practical applications the multi-stage nature should not pose a significant running time burden because most operations in each stage are analytical rather than numerical/iterative. Indeed, this computational efficiency compared with numerical or iterative approaches for fitting GLDMs, like EM, is a major advantage of subspace learning approaches such as PLDSID and our PGLDM. Finally, in terms of scalability within subspace learning approaches, PGLDM and PLDSID would scale similarly in terms of computational cost. We provide a more comprehensive explanation of this in supplementary section A.9, but briefly, we expect our method’s learning runtime to scale linearly as a function of training sample size, observation dimension, and horizon – similar to PLDSID.\n\n**References**:\n\n[1] Parima Ahmadipour, Yuxiao Yang, Edward F. Chang, and Maryam M. Shanechi. Adaptive tracking of human ECoG network dynamics.", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "ItsI9Ovbw0", "reviewer": "Reviewer_3CeQ", "summary": "In this paper, the authors propose spectral learning method for learning shared latent subspaces between multiple observations. The resulting model leverages the novel construction of Hankel matrix between observations from different processes. To promote discovery of shared subspaces, the authors propose the prioritised-extension of generliased linear dynamical modeling, which prioritises explaining the primary process using the shared subspace before learning the private subspaces specific to each process (again, priority is given to the primary process). The resulting model is demonstrated on a number of neural data benchmarks.", "questions": "See Weaknesses.", "limitations": "See Weaknesses.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is clearly written, with clear pointer to mathematical details. \n- Empirical studies is comprehensive and shows convincing results.", "weaknesses": "- Only spectral method has been proposed for learning in the proposed multi-observation GLDM. However, one unique strength of spectral methods is that the learned dynamics given spectral methods can be used as initial parameters for EM-based methods. It would be great to see if such extension is going to further improve the capability of the model.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for finding our manuscript clearly written with comprehensive and convincing experimental results. The reviewer raises a great question regarding initializing EM with our learned parameters. To the best of our knowledge, no existing EM algorithm can perform dissociation of shared vs. private dynamics and therefore there are no suitable EM algorithms that we could initialize with our parameters and use out-of-the-box.  \n\nSpecifically, there are two key features in our PGLDM learning approach: 1) we learn a model with the specific block structure defined in equation (7)  (i.e., $A_{11}$, $A_{21}$, $A_{22}$, and $A_{33}$) that allows shared and private dynamics to be dissociated within the model as distinct states $x^{(1)}$, $x^{(2)}$, and $x^{(3)}$, and 2) we learn this model in multiple consecutive stages such that the learning of shared dynamics can be prioritized in the first stage. Because existing EM methods do not allow for either the block-structure in eq. (7) or the prioritized learning of shared dynamics, we hypothesized that even if we initialize these methods with our parameters, they would not maintain the desired prioritized structure of the dynamical matrix $A$ (eqs. (6) and (7)).\n\nTo test this hypothesis we used one of our simulations from Fig 1 to learn model parameters with our method at the true shared latent dimension (i.e., $n_x=n_1$) and initialize Laplace-EM (one of our baselines) with the learned $A$ and $C_r$ parameters. We ran EM for 50 iterations and computed the eigenvalue error of the final identified modes. We compared the resulting error with the error associated with the $A$ our method had identified and the error after running Laplace-EM with random initialization. We did this experiment 10 times with different folds of 40,000 training samples, training only on the primary observation time-series. Results are presented in the table below:\n\n|Method ($n_x=n_1$)|Normalized eigenvalue error mean &plusmn; STE (log10)|\n|-|-|\n|PGLDM (stage 1)|**-2.408 &plusmn; 0.0582**|\n|Lap-EM (random init.)|-0.9351 &plusmn; 0.0074|\n|Lap-EM (PGLDM init.)|-1.0476 &plusmn; 0.0333|\n\nWe can see that initializing with our model parameters allowed EM to identify the modes more accurately than when randomly initialized, but it still had higher error than our method – as we hypothesized. There could be multiple reasons for this. Standard EM learns model parameters by maximizing the observation log-likelihood. For example, the final learned $A$ may deviate from our identified sublocks because it may be beneficial for EM to mainly learn private dynamics if they are more dominant than shared dynamics in the observation time-series. This is because standard EM has no way of prioritizing the learning of shared dynamics, i.e., our method’s second feature noted above. Also, even when EM is learning both shared and private dynamics, it will have no requirement/way to dissociate them into separate states because it lacks a block-structured approach (our method’s first feature noted above). So again EM can deviate from the clear dissociation that the sublocks identified by our method allows and instead learn states that reflect a mixture of both private and shared dynamics. Taken together, only an EM algorithm that can dissociate shared vs. private dynamics during optimization would be beneficial for being initialized with our PGLDM’s learned parameters, as also suggested by our new analysis. However, to the best of our knowledge, such an EM algorithm does not exist and the derivation/implementation of such an algorithm is better suited for future work (i.e., deriving an EM algorithm for the block structure proposed in eqs (6) and (7) or modifying the optimization cost function to prioritize shared dynamics by being multi-staged, for example).", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "5mgCZuvJDN", "reviewer": "Reviewer_m47S", "summary": "In their study \"Spectral Learning of Shared Dynamics Between Generalized-Linear Processes\", the authors introduce a multi-stage spectral algorithm for learning the parameters of a model for two generalized-linear time-series with shared latent dynamics. Assuming the latent dynamics to contain both shared and private latent dimensions, the authors show that the shared dynamics can be identified from a decomposition of a Hankel matrix formed from the time-lagged covariances between the time series. Additional (second and third stage) analyses can be used to identify the latent spaces private to either of the two time-series. Thanks to moment-conversion schemes developed here and in previous work, their spectral method can be applied to time-series with Gaussian-, Bernoulli- and Poisson-distributed observations. On synthetic data and two datasets with primate neural recordings, the authors show that their method is able to accurately identify the private and shared dynamics. The method appears to perform particularly well on identifying the shared dynamics, leading to better predictions of one time series from the other than tested alternative methods. On the neural datasets, this allows better prediction of behavioural measurements and activity in other brain areas from spiking data.", "questions": "I feel the first experiment (comparisons against Laplace-EM / PLDSID / bestLDS) could be a little better explained - are the comparisons trained only on the data of the primary time-series? The text only describes them as \"models of the primary time-series as learned with either Laplace-EM or [...]\". Is PSID trained on both primary and secondary time-series? A citation for Laplace-EM would be helpful. Also, what is the number of data points generated for the results in table 1? \n\nI am a little surprised about the bad performance of EM in Fig 1 compared to PLDSID - commonly, one would initialize a Poisson-EM fit with the much cheaper PLDSID fit and get at least a mild improvement. Are the shown results related to the 'Laplace'-part of the Laplace-EM algorithm? \n\nCould it be helpful to define $H$ as $H_r$ ? This notation would deviate from the classical literature only dealing with a single time series, but be consistent with the eventual definition of $H_z$ for step 3. It feels more natural to define $H$ as the Hankel matrix of the jointly-Gaussian stacked $[r, z]$, with $H_r$, $H_z$ and $H_{zr}$ collecting specific blocks.\n\nDid the matrices $Q$ described in line 998 conform with assumption A.1 ? \n\nIn A.7.1, why describe the norm as Froebenius rather than Euclidean when it's about vectors? \n\nChoosing red and green as important colors for the figures is not ideal for colorblind people.", "limitations": "The authors discuss limitations of their work in a dedicated section.\n\nThe authors state that by the fundamental nature of their work, they see no direct positive or negative societal consequences. I tend to disagree with the authors on this point, due to the possible use of works such as this one in Brain Computer Interfaces. Those hold promise to greatly benefit people with bodily limitations such as locked-in patients.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The suggested method adds nicely both to the subspace identification literature and a growing body of literature in computational neuroscience that is interested in identifying shared subspaces between multivariate signals. The insight that in a linear dynamical system $(x,r,z)$ structured as in equation (7), the shared dynamics $A_{11}$ are identifiable from the time-lagged cross-covariances between $r$ and $z$ is really nice.\n\nThe core derivations of the suggested method appear solid and the experiments convincingly show that the method works (which I think for most parts is all it needs to do, as for generalized-linear time-series, to my knowledge there is little direct competition). \n\nThe manuscript is written both enjoyable and succinct. In particular the methods section -- including corresponding supplementary information -- is clear and well-written (minor comments on clarity of the results section below).", "weaknesses": "As the authors remark themselves, the research question addressed here is primarily interesting for neuroscientific applications, and may hence appear somewhat niche for a larger machine learning conference. This is not a strong reservation though, in particular given the roots of the NeurIPS conferences. The contribution to the subspace identification literature should make it worthwhile to a larger audience.\n\nI am not sure if I can agree with the WLOG statement in Def 3.1 and section A.2 as it is currently written. The generality claim appears to rest on the possibility of $n_2 = n_3 = 0$ and hence $n_x = n_1$ - but in this case, what is the point of this work? I feel that the authors first should introduce the structure of shared and private latent spaces (which is not yet established by the preceding section 3.1) with the respective dimensionality split $n_x = n_1 + n_2 + n_3$. Then, given that assumption of the structured latent space with fixed $(n_1, n_2, n_3)$, equation (7) would be general.\nThe assumption A.1 is important for the algorithm and should be stated in the main text, ideally in definition 3.1 or in a small remark following it.\n\nThe multi-stage nature of the procedure feels somewhat prone to accumulation of errors -- since this is a spectral method this will hold particularly for small sample sizes. More specifically, I am a somewhat worried about numerical stability following the subtractions and subsequent matrix decompositions in equations (10) and (22) - any comments on this in the presence of estimation noise, in particular if one of the two terms shown in eq (19) to make up H dominates the other, e.g. in terms of Froebenius norm?\n\nThere is a typo in lines 822 (dropped a subscript '11' from one of the occurrences of the first block of A), and a missing subscript 'r' on matrix $C_r$ in line 967.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their close read of our manuscript and for finding it enjoyable with a well-written methods section. Below we address comments/questions inline.\n\n> Feedback on Def 3.1, section A.2, & assumption A.1\n\nWe thank the reviewer for their suggestions. We clarify that the WLOG in Def 3.1 and section A.2 was to suggest that our algorithm could cover the learning of a GLDM for the single observation time-series as a special case (i.e., using stage 2 only with n1=n3=0, and n2=nx). We now realize this was unclear and have removed the WLOG statement in section 3.1 and revised section A.2 to make our intended meaning clearer. We will also move Assumption A.1 to section 3.1.\n\n> accumulation of errors in the presence of estimation noise\n\nWe thank the reviewer for their great question. There are two aspects of the learning we wish to comment on.\n\nFirst there is a signal-to-noise (SNR) consideration that is fundamental to learning methods in general. In stage 1 our method requires high SNR in the cross-correlations between future secondary time-series observations and past primary time-series observations. If most of the signal present in the primary time-series is not attributable to the shared latent states (i.e., the residual Hankel in eq (19) dominates), then most of $H_{zr}$’s singular values will be small and stage 1 may have greater estimation error. If the reverse situation holds and the shared latent states explain most of the signal in the primary time-series, then the residual Hankel $H_{r}^{(2)}$ will mostly have small singular values, possibly resulting in estimation errors during stage 2. Inspection of the singular values prior to model parameter extraction can, however, help guide method usage. For example, if the singular values of $H_{zr}$ are small, then this may indicate that the two time-series don’t have shared dynamics and only stage 2 of the method is needed (i.e., standard GLDM). If the singular values of the residual Hankel are small, then this may indicate that almost all dynamics are shared and only stage 1 is needed. \n\nA second consideration is a numerical one that could, as the reviewer stated, result in an accumulation of errors in downstream stages. If during stage 1 the modes are identified inaccurately because there is minimal shared dynamics between the two time-series and/or the training sample size is too small, then this estimation error could impact the computation of the residual matrix $H_{r}^{(2)}$, introducing error in stage 2. Examining the singular values of $H_{zr}$ may, however, also help avoid this situation. Although the multi-stage learning can lead to error accumulation in some situations, it comes with the benefit that our method has the ability to more accurately identify the shared dynamics (when they exist) in stage 1, compared to existing GLDM methods. Finally, Figs 1b, 2b, and 3b suggest that the impact of error accumulation is not severe and can be situation-dependent; for example, self-prediction performance of PGLDM reaches that of PLDSID in 1b, reaches very close to it in 2b, and exceeds it in 3b. We will add a supplementary section discussing these points and add error accumulation to the Limitations section as a potential downside.\n\n> clarification on the first experiment\n\nThank you for the feedback. We’ve revised the manuscript (sections 4.1 and A.6.1) to include more details regarding the first experiment. PSID and our method were trained on both the primary and secondary time-series because these algorithms can consider both time-series for training. All other baselines used only the primary time-series, as these methods only model a single data source. We’ve now included a citation in Table 1 to the GitHub library for the Laplace-EM algorithm that we used (previously cited only as a footnote). Lastly, we simulated 25600 timesteps for each configuration.\n\n> bad performance of EM in Fig 1\n\nWe thank the reviewer for their great question. It is possible that the global Laplacian approximation underlying the EM algorithm in [1] may explain its performance in Fig 1. In this widely-used EM implementation [1], the global posterior (i.e., the probability of the latents $x$ conditioned on all observations $y_{1:T}$) is approximated as a Gaussian distribution so that it can be computed. It may be the case that the accuracy of the Gaussian approximation is poor for the posterior related to the Poisson observations in Fig 1, impacting overall model performance.\n\n> Structure of synthetic $Q$ matrices\n\nWe thank the reviewer for the great question. For most of our simulations the $Q$ matrices did conform to assumption A.1. We deviated from the assumption only for the Poisson-Poisson simulations used in Table 1. However, as shown by the results in Table 1, this did not significantly impact our method's identification of shared dynamical modes. In most cases we simulated the block $Q$ structure by simulating the dynamics private to the secondary observation time-series as a separate latent dynamical model. Because latent $x_{3}$ is completely decoupled from $x_{1}$ and $x_{2}$ in our model definition (eqs (7) and (13), assumptions A.1 and A.2), we chose to generate one dynamical model corresponding to the latents $x_{1}$ and $x_{2}$ and a separate dynamical model corresponding to $x_{3}$. This approach generates a block $Q$ but only works for Gaussian observations, which is why the Poisson-Poisson case was different. We will add this clarification to section A.6.1.\n\n> Additional comments\n\nWe have fixed the typos in lines 822 & 967, renamed $H$ to $H_r$, renamed the Frobenius norm to l2-norm, and will plan to adjust our plot colors in the final manuscript appropriately. We thank the reviewer for their recommendations. We also thank the reviewer for their assessment of our work and its applicability in the BCI domain. We will include this point in the manuscript as a potential benefit.\n\n**References**:\n\n[1] Linderman lab SSM library", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "asyPq3Hhlk", "reviewer": "Reviewer_95d5", "summary": "The work proposes a new linear dynamical system model that aims to model interaction between two time series by explicitly modeling their shared and private dynamics. Furthermore, authors incorporate generalized linear models in their observation model of the proposed dynamical system to handle different kinds of time series typically encountered in neuroscience (e.g., count-time series, real-valued time series etc.). Authors then extend standard covariance-based subspace system identification to learn the parameters of the proposed model and compared the performance of their proposed model with existing baselines on two experimental datasets and a simulation study.", "questions": "- The structure of the proposed $A$ matrix is a bit unintuitive and does not seem to satisfy intuitive properties one would expect from a dynamical system modeling shared and private latents. Intuitively, I would have expected the structure of the $A$-matrix to have the forms of these kinds: $A=\\begin{bmatrix}A_{11}&0&0\\\\\\\\ 0&A_{22}&0 \\\\\\\\ 0&0&A_{33}\\end{bmatrix}$, where the dynamics of the shared latent $x^{(1)}$  and the dynamics of the private latents $x^{2}$ and $x^{3}$ do not affect each other or the structure $A=\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\\\\\ A_{12}&A_{22}&0 \\\\\\\\ A_{31}&0&A_{33}\\end{bmatrix}$, where the dynamics of the private latents $x^{3}$ and $x^{2}$ would only affect each other through the shared latent $x^{(1)}$. The structure proposed by the authors, i.e. $A=\\begin{bmatrix}A_{11}&0&0\\\\\\\\ A_{12}&A_{22}&0\\\\\\\\ 0&0&A_{33}\\end{bmatrix}$ has an inherent asymmetry in-built to its dynamics where the dynamics of the second private latent $x^{(2)}$ is affected by the shared latent state but the private latent $x^{(3)}$'s dynamics are only affected by itself. Why is this modeling choice justified? Why dynamics of only one of the private latent can be affected by the shared latent. Is there a domain-specific reason?   \n- Why does PGLDM always perform worse at self-prediction compared to other baselines (Figure 1 (b), 2(b), 3(b))?\n- I am also curious to know that how the results look when decoding the activity V2 from V1, and V1 self-prediction using the system estimated by PGLDM in Section 4.3", "limitations": "Authors have provided a section discussing some limitations. A major limitation of this work that authors do not mention is the comparison with non-linear models which are also being used for decoding neural activity (such as neural data transformer by Ye & Pandarinath'21). I do not think a experimental comparison between these transformer architecture and the linear models proposed in this work is crucial, as the goal of these models have widely diverged with linear models being used more for understanding underlying neural dynamics, whereas the transform models being used for boosting decoding performance. Regardless, I think a mention of these non-linear methods should be made in the limitation sections along with the acknowledgement of a lack of comparison with these methods.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is very well-written and the ideas are introduced in a natural progression making them easy to follow.\n- Authors identify a gap in literature where existing methods do not model multiple time-series while explicitly accounting for shared and private dynamics when the observed time series are not real-valued.\n- Authors propose an interesting dynamical system for modeling two different time series containing shared and private latents, whose parameters can be analytically estimated by extending existing covariance-based subspace system identification algorithms. \n- The proofs provided in the paper are generally correct.\n- Authors show improvement in the models' ability to predict one time series properties from the other time series when using their proposed dynamical system model which models the dynamics of the two time series jointly.", "weaknesses": "- The proof of the claim that the any general $A$, $C_r$, and $C_z$ can be decomposed into the block form structure described in equation (7) is incomplete. The proof provided in Appendix A.2 uses a well-known theorem, which states that for a non-fully observable linear dynamical system characterized by the system matrices $(A,C)$, there exists a non-singular linear transformation $T$ such that the new linearly transformed system having the new system matrices $A'=T^{-1}AT$ and $C'=CT$ have the following block structure:\n$ \nTA^{-1}T = \\begin{bmatrix} A_{11} & A_{12} \\\\\\\\ 0 & A_{22} \\end{bmatrix}, CT = \\begin{bmatrix}0& C_2\\end{bmatrix}\n$, where the pair $(A_{22}, C_2)$ are observable. Authors apply this theorem to the **special case** when the dimension of $x^{(3)}$ is zero, to reduce the matrices $A, C_z$ and $C_r$ to the desired block diagonal form, which is correct. But then, authors state that a similar argument works when the dimension of $x^{(3)}$ is non-zero. I am not sure how this proof of the special case of $x^{(3)}$ having zero-dimension can be extended to the general case of $x^{(3)}$ having non-zero dimension. The argument that the general matrices $A=\\begin{bmatrix}A_{11}& A_{12}& A_{13}\\\\\\\\ A_{21}&A_{22}& A_{23}\\\\\\\\ A_{31}&A_{32}&A_{33} \\end{bmatrix}$, $C_z = \\begin{bmatrix} C_z^{(1)}&C_z^{(2)}&C_z^{(3)}\\end{bmatrix}$, and $C_r=\\begin{bmatrix}C_r^{(1)}&C_r^{(2)}&C_r^{(3)} \\end{bmatrix}$ can always be reduced to the block-diagonal shown in equation (7) does not fall out as a trivial consequence of the theorem stated above, and requires an explicit proof. Naively applying the above theorem using the dimension of the observability matrix pair $(A, C_z)$ would show an existence of $T$ such that the new transformed system has the matrix with the following block structure: $A=\\begin{bmatrix}A_{11}& A_{12}& A_{13}\\\\\\\\ A_{21}&A_{22}& A_{23}\\\\\\\\ 0&0&A_{33} \\end{bmatrix}$, $C_z = \\begin{bmatrix} 0&0&C_z^{(3)}\\end{bmatrix}$.  I would ask authors to provide a rigorous mathematical proof showing that there always exist a non-invertible matrix $T$ such that a general $A, C_z$ and $C_r$ can be reduced to the block-diagonal form shown in equation (7). In the event, authors are unable to provide a proof, I would ask authors to modify the manuscript accordingly to reflect that the block diagonal structure is an assumption on the system model and remove the phrases \"without loss of generality\". In this review, I am going to continue on the assumption that the general matrices $A, C_r$ and $C_z$ cannot be reduced to the block-diagonal form (until shown otherwise by the authors through a rigorous mathematical proof) and the block-diagonal form is an assumption on the dynamical systems the authors are trying to estimate.  \n- The claims of the papers that the proposed model can be learnt for any generalized-linear process are overstated. To be able to learn the proposed models' parameters (described in equation (6)), the moment conversion equations described in section 3.2.3 are needed. Authors only provide the moment conversion equations for Gaussian, Poisson, and Bernoulli distributions. It is unclear how the model parameters can be learnt when, for example, $\\mathcal{P}_{y|r}(y_k; g(r_k))$ is an exponential, gamma, or inverse-gaussian distribution. It may be the case that a general moment conversion equation may not exist for any arbitrary exponential distribution. Accordingly, I would ask the authors to town down their claims or provide a general moment conversion equation for any arbitrary exponential family distribution.\n\nMinor\n\n- Please do not use pythonic notation (without defining it first) for mathematical vector operations. For example in equation (11), R.H.S. of the equation is dividing a matrix by a vector. That is not a well-defined operation.  Similarly, equation (12) the operation $\\mu_{t_{f_m}}\\mu_{y_{p_n}}$ does not make sense, it should be $\\mu_{t_{f_m}}\\mu_{y_{p_n}}^T$. Please also state that $\\ln(\\cdot)$ in those operations is applied elementwise, as there exists well defined ways to define matrix functions of $\\ln(\\cdot)$ too.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for finding that the manuscript is well-written and addresses a gap in the literature. We also thank them for providing very comprehensive and helpful feedback. Below we reply to outstanding questions/concerns inline.\n\n> Generality claim on block form structure described in eq (7)\n\nWe thank the reviewer for their feedback. The reviewer is exactly correct in all their points and we have revised the manuscript and Appendix A.2 to reflect these points. Briefly, as the reviewer points out, our proof on model generality provided in Appendix A.2 only applies to the case of $n_3=0$ (when learning the private dynamics of the secondary time-series is not of interest). When extending our method to the case of $n_3 > 0$ (i.e., optional stage 3), we add the following assumptions that are all motivated by having well-defined *private* dynamics but are not necessarily general:\n1) We assume that the private dynamics of z ($x^{(3)}$) do not drive the shared dynamics of y and z ($x^{(1)}$), so as not to leak private z dynamics into shared dynamics. Thus we take $A_{13}=0$.\n2) We assume that the private dynamics of the two signals ($x^{(2)}$ and $x^{(3)}$) do not drive each other, and for this reason we take $A_{23}=0$ and $A_{32}=0$. \n3) Our assumption that $A_{31}=0$ was motivated by the definition of private dynamics in the context of our application. In our use case, the primary time-series ($y_k$) is used as the only observation during inference after the model is trained with both primary and secondary ($z_k$) time-series. In this setup, if $A_{31}$ was nonzero, the primary time series would be able to infer/predict, through $A_{31}$, private latents $x^{(3)}$ that are predictive of the private secondary time-series dynamics, undermining the notion of privacy. We do not encounter this situation with a non-zero $A_{21}$ and private primary dynamics $x^{(2)}$ because we do not use the secondary time-series as observation to estimate the primary time-series. \n\nWe have revised Appendix A.2 and the main text to clarify these additional assumptions in the case of $n_3>0$, their motivation in terms of having well-defined “private” latents, and the fact that the complete generality in the block formulation of A is only for the $n_3=0$ case (stages 1 and 2). Future work interested in symmetrical inference applications may find ways to forgo some of the above assumptions. \n\n> Structure of the proposed $A$ matrix\n\nWe thank the reviewer for their great question. We clarify a few points:\n- Our approach is also amenable to the reviewer’s first proposed structure $A=\\begin{bmatrix}A_{11}&0&0\\\\\\\\0&A_{22}&0\\\\\\\\0&0&A_{33}\\end{bmatrix}$. This is because in stage 2 we can extract $A_{22}$ directly by computing the least-squares solution $\\underline{\\Delta}^{(2)} \\overline{\\Delta}^{(2)\\dagger}$ (section 3.2.2 and eq. (24)). However, we chose to keep $A_{21}$ to maintain the most general form (which in the case of $n_3=0$ is well-supported by theorem 3.8 [1], as discussed previously).\n- The second proposed structure $A=\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\\\\\A_{21}&A_{22}&0\\\\\\\\A_{31}&0&A_{33}\\end{bmatrix}$ would not be suitable for our use case because unshared/private latents $x^{(2)}$ and $x^{(3)}$ would contribute to future timesteps of shared latents $x^{(1)}$ (via $A_{12}$ and $A_{13}$) – coupling shared and unshared; if we allow private dynamics to contribute to the future of shared dynamics, then this means they are not fully private. In our formulation we chose to keep the shared states completely isolated in their recurrent dynamics.\n- Arguably the most general form would have been to allow a nonzero $A_{31}$ in our formulation:  $A=\\begin{bmatrix}A_{11}&0&0\\\\\\\\A_{21}&A_{22}&0\\\\\\\\A_{31}&0&A_{33}\\end{bmatrix}$. However, as noted above, we chose to assume $A_{31}$ = 0 because otherwise the primary time-series would be informative/predictive of private dynamics in the secondary time-series, undermining the notion of privacy. \n\n> PGLDM self-prediction performance\n\nThank you for the question. When using stage 1 only, our algorithm dedicates all the latent state dimensions ($n_1$) to the prioritized identification of shared dynamics. By doing so, the identified models are able to more accurately learn the shared dynamics and more accurately decode the secondary time-series. However, this means that self-prediction in stage 1 may suffer because the dominant dynamics in the primary time-series may be private rather than shared. After adding stage 2, however, PGLDM starts to learn the private dynamics and so its self-prediction performance improves in Figs 1b, 2b, and 3b, reaching very close to self-prediction of PLDSID once given enough latent state dimensions.\n\n> V2 decoding from V1; V1 self-prediction\n\nWe thank the reviewer for their question. We have added a figure of the results in the PDF attached to our general response, which we will also include in the appendix. We originally performed our analysis only in the V2-V1 direction because prior work has shown that the feedback direction, that is the V2 (past) to V1 (future) direction, generally exhibited higher correlations between the two time-series as compared to the feedforward direction [2]. This might also explain the differences that we see in the results.\n\n> Additional comments\n\nWe will include a discussion regarding nonlinear methods in the Limitations section as suggested. We have corrected the use of Pythonic notation indicated by the reviewer. We have also revised the text in lines 297-303 to tone down the claims about applicability to generalized-linear processes and now emphasize the need for a computationally tractable moment conversion equation. We thank the reviewer for the feedback.\n\n**References**:\n\n[1] Tohru Katayama. Subspace Methods for System Identification.\n\n[2] Semedo, J.D., Jasper, A.I., Zandvakili, A. et al. Feedforward and feedback interactions between visual cortical areas use different population activity patterns.", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for taking the time to review our submission and for providing helpful feedback, suggestions, and discussion regarding our work. We were encouraged to hear that reviewers found our manuscript “enjoyable” (reviewer m47S) and “well-written” (reviewers Hfka, m47S, and 95d5), with a “comprehensive” (reviewer 3CeQ) experimental section demonstrating that our “interesting” (reviewer 95d5) and “novel” (reviewer Hfka) approach works.
 
We provide responses to each reviewer’s comments and questions inline in each of our rebuttals below. We have also made revisions to the manuscript as needed to further address reviewer comments. Lastly, we performed a few new analyses for the manuscript to address some reviewer comments, and are including these as a PDF with a revised table 2 (with new running time analysis results), a revised figure 3 (additional barplots for V2-V1 analysis), and a new supplementary figure to complement figure 3 (V1-V2 analysis). Tables/figures are presented in that listed order. The color scheme in the attached figure, as well as all figures in the manuscript, will be updated in the final version to improve accessibility for color blind individuals, as advised by reviewer m47S.

Meta-review: This paper describes a subspace ID method for characterizing generalized dynamical models.  All four reviewers felt that it was above the threshold for acceptance, and I'm pleased to report that it has been accepted to NeurIPS.  Congratulations!  Please revise the manuscript according to the reviewer comments and discussion points.
Decision: Accept (poster)