Spectral Learning of Shared Dynamics Between Generalized-Linear Processes

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 20:13:04

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.

Main paper: 
 introduction: 
 Modeling the shared dynamics between temporally-structured observations with different statistical properties is useful across multiple application domains, including neuroscience and biomedical engineering (D'mello & Kory, 2015; Lu et al., 2021) . However, building models of the dynamic relation between such signals is challenging for two key reasons. First, continuous-and discretevalued observations exhibit different statistics, which the modeling approach must appropriately reconcile. Second, residual (i.e., unshared or unique) dynamics present in each observation stream can obscure and confound modeling of their shared dynamics (Allen et al., 2019; Stringer et al., 2019; Sani et al., 2021) . Thus, the modeling approach also needs a way to accurately dissociate and prioritize identification of the shared dynamics. Current analytical methods do not simultaneously enable both of these capabilities, which is what we address here.Linear dynamical state-space models (SSMs) are a commonly used framework for modeling dynamics using a low-dimensional latent variable that evolves over time (Paninski et al., 2010; Macke et al., 2015; Newman et al., 2023) . Even though the past decade has seen an increased use of artificial neural networks and deep learning methods for training dynamical models of time-series data (Pandarinath et al., 2018; Hurwitz et al., 2021; Kramer et al., 2022; Schneider et al., 2023) , analytical SSMs still remain widely popular due to their interpretability and broad applicability both in scientific investigations and in real-time engineering applications (Kao et al., 2015; Aghagolzadeh & Truccolo, 2016; Lu et al., 2021; Yang et al., 2021; Newman et al., 2023) . For Gauss-Markov 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 and dissociation of intrinsic vs input-driven activity (Van Overschee & De Moor, 1996; Katayama, 2005; Sani et al., 2021; Galgali et al., 2023; Vahidi et al., 2023) . These methods, however, are not applicable to generalized-linear processes with non-Gaussian observations. While there has been work extending SSID to generalized-linear processes, such as Poisson and Bernoulli observations (Buesing et al., 2012; Stone et al., 2023) , these methods only learn the dynamics of a single observation time-series rather than modeling shared dynamics between two time-series (see section 2.1). Finally, prior multimodal learning algorithms do not explicitly tease apart the shared vs. residual (disjoint) dynamics in a predictor (primary) time-series, but instead model the collective dynamics of two modalities in the same latent states (Abbaspourazad et al., 2021; Kramer et al., 2022; Ahmadipour et al., 2023) .Here we fill these methodological gaps by deriving a novel covariance-based SSID learning algorithm that (1) is applicable to generalized-linear processes, and ( 2) is capable, with its two-staged learning approach, of identifying with priority the shared dynamics between two processes before modeling residual (predictor-only) dynamics. To illustrate the method, we focus on the specific case of modeling Poisson-distributed discrete time-series while dissociating their shared dynamics with Gaussian-distributed continuous observations, which is of particular interest in neuroscience. However, we emphasize that our method can be extended to other output distributions in the generalizedlinear model family (section 5). We show that our method successfully dissociated the shared dynamics between Poisson and Gaussian observations both in simulations and on a public non-human primate (NHP) dataset of discrete population spiking activity recorded during continuous arm movements (O'Doherty et al., 2017) . Further, compared with existing Poisson SSID methods, our method more accurately decoded movements from Poisson spiking activity using lower-dimensional latent states. Lastly, we discuss limitations and propose potential solutions and future research directions. 
background: 
 Our method provides the new capability to dynamically model Poisson observations, while prioritizing identification of dynamics shared with Gaussian observations. We first review the existing SSID method for modeling Poisson observations, which serves as our baseline, as well as standard covariance-based SSID, to help with the exposition of our method in section 3.A Poisson linear dynamical system (PLDS) model is defined asx k+1 = Ax k + w k r k = C r x k + b y k | r k ∼ Poisson(exp(r k )) (1)where x k ∈ R nx is the latent state variable and y k ∈ R ny corresponds to discrete (e.g., neural spiking) observations which, conditioned on the latent process r k , is Poisson-distributed with a rate equal to the exponential of r k (i.e., log-rate). Finally, N (w k ; 0, Q) is state noise and b is a constant baseline log-rate. The PLDS model is commonly used for modeling Poisson process events, such as neural spiking activity (Smith & Brown, 2003; Truccolo et al., 2005; Lawhern et al., 2010; Buesing et al., 2012; Macke et al., 2015) . Buesing et al. (2012) developed an SSID algorithm, termed PLDSID, to learn the PLDS model parameters Θ = (A, C r , b, Q) given training samples y k and hyperparameter n x corresponding to the latent state dimensionality.There exist standard covariance-based SSID algorithms (section 2.2) that can learn the parameters of a latent dynamical system given a future-past Hankel matrix, H, constructed from the crosscovariances of the system's linear observations as (Van Overschee & De Moor, 1996; Katayama, 2005 )EQUATIONwhere the integer i denotes the user-specified maximum temporal lag (i.e., horizon) used to construct H and Λ rτ := Cov(r k+τ , r k ) is the τ -th lag cross-covariance for any timepoint k, under time-stationary assumptions. Such covariance-based SSID algorithms, however, are not directly applicable to Poisson-distributed observations (section 2.3). This is because the log-rates r k that are linearly related to the latent states in equation ( 1) are not observable in practice -rather, only a stochastic Poisson emission from them (i.e., y k ) is observed. As a result, the second moments constituting H (i.e., Λ rτ ) cannot be directly estimated. The critical insight by Buesing et al. (2012) was to leverage the log link function (i.e., exp -1 ) and the known conditional distribution y k |r k to compute the first (µ r ± ) and second (Λ r ± ) moments of the log-rate r k from the first (µ y ± ) and second (Λ y ± ) moments of the discrete observations y k . The ± denotes that moments are computed for the future-past stacked vector of observations r ± := r T f r T p T and y ± := y T f y T p T , whereµ r ± := E[r ± ] µ y ± := E[y ± ] Λ r ± := Cov(r ± , r ± ) Λ y ± := Cov(y ± , y ± ).To compute moments of the log-rate, Buesing et al. (2012) derived the following moment conversionEQUATIONwhere m ̸ = n correspond to different indices of the first and second moments of the future-past stacked observation vectors r ± and y ± , and n, m = 1, • • • , Kn y where K is the total number of time points. With the first and second moments computed in the moment conversion above, the baseline log rate b parameter is read off the first n y rows of µ r ± and the Hankel matrix, H, is constructed as per equation ( 2). From here, it is possible to proceed with the standard covariancebased SSID algorithm for Gauss-Markov models using H, as outlined next.Given an H matrix, covariance-based SSID first decomposes H into a product of observability (Γ r ) and controllability (∆) matrices as (Van Overschee & De Moor, 1996; Katayama, 2005 )H SVD = Γ r ∆ =     C r C r A . . . C r A i-1     A i-1 G • • • AG G (4)where G := Cov(x k+1 , r k ). The factorization of H is done by computing a SVD of H and keeping the top n x singular values and corresponding singular vectors. Note that the rank of H must be at least n x in order to identify a model with a latent dimension of n x . Thus, the user-specified horizon i must satisfy i × n y ≥ n x . From the factors of H, C r is read off as the first n y rows of Γ y and A is learned by solving Γ r = Γ r A, where Γ r and Γ r denote Γ r from which the top or bottom n y rows have been removed, respectively. This optimization problem has the following closed-form leastsquares solution A = Γ † r Γ r , with † denoting the pseudo-inverse operation. Discussion regarding learning the state noise covariance model parameter Q is postponed to section 3.2.3 below.At a high-level, there exist two common approaches for subspace identification (Van Overschee & De Moor, 1996; Katayama, 2005) : (i) covariance-based methods (e.g., Buesing et al. (2012) ; Ahmadipour et al. (2023) ) that aim to learn all model parameters based on the second-order statistics of the observations and not the observation time-series directly, and (ii) projection-based methods (e.g., Sani et al. (2021) ; Vahidi et al. (2023) ) that make direct use of the observation time-series via linear projections. Projection-based methods are often used to model Gaussian time-series but are not applicable to Poisson observations, which instead require a covariance-based approach since the latent log-firing rates (r in equation ( 1)) are unobserved. To achieve our aim (i.e., modeling a generalized-linear process with prioritized identification of shared dynamics), we need to develop a novel covariance-based subspace identification algorithm, which presents the following challenges:1. Covariance-based methods, including PLDSID (Buesing et al., 2012) , do not guarantee a valid set of parameters that satisfy the positive semidefinite covariance sequence requirements (Van Overschee & De Moor, 1996) . To address this challenge, we use the optimization approach outlined in section 3.2.3 to ensure validity of noise statistics and enable inference from the learned model.2. We could not rely on time-series projections to isolate the residual predictor process dynamics at the beginning of the second stage (Sani et al., 2021) . As a result, we derived a new least-squares problem for learning the components of the state transition matrix A corresponding to the unique predictor process dynamics (section 3.2.2), without changing the shared components learned in the first stage (section 3.2.1). By doing so, prioritized learning of shared dynamics is preserved. 
modeling shared dynamics between poisson and gaussian observations: 
 The PLDS model (equation ( 1)) is for modeling Poisson observations on their own rather than with Gaussian observations. To enable this capability, we write the following Poisson-Gaussian linear dynamical system model     x k+1 = Ax k + w k z k = C z x k + ϵ k r k = C r x k + b y k | r k ∼ Poisson(exp(r k )) (5)where z k ∈ R nz represents continuous observations (e.g., arm movements), ϵ k represents their noise (either white, i.e., zero-mean temporally uncorrelated Gaussian noise, or colored, i.e., zeromean temporally correlated Gaussian noise), and y k ∈ R ny represents the discrete observations (i.e., neural spiking). Further, we introduce a block structure to the system (Sani et al., 2021) that allows us to dissociate shared latents from those that drive the Poisson observations only. Specifically,A = A 11 0 A 21 A 22 C z = C (1) z 0 C r = C (1) r C (2) r x = x (1) x (2) (6)where x(1) k ∈ R n1 corresponds to latent states that drive both z k and y k , and x(2) k ∈ R nx-n1 to states that only drive y k . The parameter G can also be written in block partition format such thatG = E x (1) k+1 x (2) k+1 r T k -E x (1) k+1 x (2) k+1 E[r k ] T = E[x (1) k+1 r T k ] E[x (2) k+1 r T k ] - E[x (1) k+1 ]E[r k ] T E[x (2) k+1 ]E[r k ] T = G (1) G (2) .Our method, termed PG-LDS-ID (Poisson-Gaussian linear dynamical system identification), learns the model parameters, i.e., Θ ′ = (A, C z , C r , b, Q), given training samples y k and z k and hyperparameters n 1 and n 2 = n x -n 1 denoting the shared and residual latent dimensionalities. Selection of appropriate values for hyperparameters can be done with cross-validation (see appendix A.3).PG-LDS-ID uses a two-staged learning approach to model Poisson time-series while prioritizing identification of the dynamics shared with Gaussian observations. During stage 1, shared dynamics are learned using both observations. In stage 2, any residual dynamics in the predictor observations 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 predictor dynamics only if there are enough latent states to explain the shared dynamics (full derivation in appendix A.1). Note, predictor refers to the data stream whose modeling is of primary interest (the Poisson observations here) and that is used to predict the secondary data stream (the Gaussian observations). For example, Poisson observations are the predictor within the context of decoding continuous behaviors from discrete population spiking activity. The roles can be swapped without loss of generality and the designation is made clear in equation ( 8) below.In the first stage, our algorithm identifies the parameter set corresponding to the shared dynamical subspace, (A 11 , C(1) r , C z , b), given hyperparameter n 1 and using both Gaussian and Poisson observations, z k and y k . To do this, we first compute a moment conversion to estimate joint moments of z k and r k from the joint moments of the observed signals z k and y k with the following equation derived using the conditional statistical properties of both observations (see appendix A.1.4)EQUATIONNext, we use these moments to construct a Hankel matrix between future continuous observations and past log-rates of the discrete observationsEQUATIONwith r p defined as in equation ( 2). Although equation ( 8) uses the same horizon for both observations, in practice we implement the method for a more general version with distinct horizon values i y for the discrete observations and i z for the continuous observations, resulting in H zr ∈ R iz * nz×iy * ny . 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.3.1). Further, and importantly, by using z as the future observations in the Hankel matrix, we learn a dynamical model wherein Poisson observations y can be used to predict the Gaussian observations z. After constructing H zr , we decompose it using SVD and keep the top n 1 singular values and their corresponding singular vectorsH zr SVD = Γ z ∆ (1) =     C z C z A 11 . . . C z A i-1 11     A i-1 11 G (1) • • • A 11 G (1) G (1) (9)where n 1 is the user-specified dimensionality of the shared latent states x(1) k , Γ z denotes the observability matrix for the continuous observations, and ∆ (1) denotes the controllability matrix associated with the shared latent states (defined as in equations ( 4) and ( 17) in the appendix). At this point, we extract C z by reading off the first n z rows of Γ z . To extract C(1) r we first form H per equation ( 2) and extract the observability matrix for r associated with the shared latent dynamics, Γ(1) y , by right multiplying H with the pseudoinverse of ∆ (1)H∆ (1) † = Γ (1) r =       C (1) r C (1) r A 11 . . . C (1) r A i-1 11       .We then read C(1) r from the first n y lines of Γ(1) r (defined as in equation ( 20) in the appendix). The baseline log rate b is read off the first n y rows of µ r ± computed in the moment conversion from equation (3). Lastly, to learn the shared dynamics summarized by the parameter A 11 , we solve the optimization problem ∆ (1) = A 11 ∆( 1) where ∆ (1) and ∆(1) denote ∆ (1) from which n y columns have been removed from the right or left, respectively. The closed-form least-squares solution for this problem is 1) ) † . This concludes the learning of the desired parametersEQUATIONEQUATIONr , C z , b), given hyperparameter n 1 , in stage 1.After learning the shared dynamics, our algorithm can learn the residual dynamics in the predictor observations that were not captured by x(1)k . Specifically, we learn the remaining parameters from equation ( 6):EQUATIONr ), with hyperparameter n 2 = n x -n 1 determining the unshared latent dimensionality. To do so, we first compute a "residual" Hankel matrix, H (2) , using Γ(1) r and ∆ (1) from stage 1 and decompose it using SVD, keeping the first n 2 singular values and vectorsH (2) = H -Γ (1) r ∆ (1) SVD = Γ (2) r ∆ (2) . (10) With C (2)r , which corresponds to the first n y rows of Γ(2) r , we construct C r = C(1) r(2) r.We then use ∆ (2) to form the controllability matrix ∆ as the concatenation of ∆ (1) and ∆ (2) (derivation in appendix A.1):∆ = A i-1 G • • • AG G = ∆ (1) ∆ (2) . Given ∆, we next extract [A 21 A 22 ] by solving the problem ∆ (2) = [A 21 A 22 ] ∆ where ∆ (2) := [A 21 A 22 ] A i-2 G • • • [A 21 A 22 ] G , ∆ := A i-2 G • • • G .A = A 11 0 A 21 A 22, we now have all model parameters,(A, C r , C z , b), given hyperparameters n 1 and n 2 , except state noise covariance Q.Standard SSID algorithms (e.g., section 2.2) learn linear SSMs of the following formx k+1 = Ax k + w k r k = C r x k + v k (11)where the new term N (v k ; 0, R) corresponds to observation noise. State noise, w k , and observation noise, v k , can have a non-zero instantaneous cross-covariance S = Cov(w k , v k ). SSID in general does not assume any restrictions on the noise statistics. However, the Poisson observation model (equations ( 1) and ( 5)) has no additive Gaussian noise for r k and instead exhibits Poisson noise in y k , when conditioned on r k . This means that v k = 0 in equation ( 5), and thus R = 0 and S = 0. Imposing these constraints is important for accurate parameter identification for Poisson observations, but was not previously addressed by Buesing et al. (2012) . Thus, we require our algorithm to find a complete parameter set Θ ′ that is close to the learned (A, C r , C z , b) from the two stages in sections 3.2.1 and 3.2.2 and imposes the noise statistic constraints R = 0 and S = 0. To do this, inspired by Ahmadipour et al. (2023) , we form and solve the following convex optimization problem to satisfy the noise statistics requirements minimizeΛx ∥S(Λ x )∥ 2 F + ∥R(Λ x )∥ 2 F such that Λ x ⪰ 0, Q(Λ x ) ⪰ 0, R(Λ x ) ⪰ 0 (12)where Λ x := Cov(x k , x k ) denotes the latent state covariance and the following covariance relationships, derived from equation ( 11) (Van Overschee & De Moor, 1996) , holdEQUATIONThis 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.3). It also enables state prediction (see appendix A.4). Combining the previously found parameters and the matrix Q that corresponds to the minimizing solution Λ x of equation ( 12), we have the full parameter setΘ ′ = (A, C r , C z , b, Q).We used Python's CVXPY package to solve the semidefinite programming problem defined in equation ( 12) (Diamond & Boyd, 2016; Agrawal et al., 2018) . For all of our comparisons against baseline, we learned the noise statistics associated with PLDSID's identified parameters using this approach, keeping the rest of the algorithm the same. 
shared dynamics are accurately identified in simulations: 
 We simulated Poisson and Gaussian observations from random models as per equation ( 5) to evaluate how well our method identified the shared dynamics between the two observations. All state and observation dimensions were randomly selected and the corresponding system parameters generated to simulate stable and slow-decaying dynamics (see appendix A.5.1). We computed two performance metrics: 1) the normalized eigenvalue error between ground truth and identified shared dynamical modes (i.e, the eigenvalues of A 11 in equation ( 6)), and 2) the predictive power of the model when using discrete Poisson observations to predict continuous Gaussian observations in a held-out test set. This second metric allowed us to test our hypothesis that PG-LDS-ID's explicit modeling of the shared subspace improved decoding of Gaussian observations from Poisson observations compared with PLDSID (Buesing et al., 2012) . To compute the first metric for PLDSID, which does not explicitly model shared dynamics, we needed to select the n 1 modes identified from the Poisson time-series only that were the most representative of the Gaussian time-series. To do so, we first trained PLDSID on Poisson observations and extracted the latent states. Then, we sorted these learned latent states based on their accuracy in predicting the Gaussian observations (appendix A.4). We computed the eigenvalues associated with the top n 1 most predictive latent states, which we considered as the shared modes identified by PLDSID. We computed the normalized eigenvalue error as For configurations wherein learned n x is smaller than true n 1 , we substituted missing modes with 0 prior to computing the normalized error.|Ψ true -Ψ id | F / |Ψ true | F ,In our first simulation experiment, we generated 50 random systems and studied the effect of training set size on learning. We used 1e2, 1e3, 1e4, 1e5 or 1e6 samples to train models and tested them on 1e6 samples of independent held-out data (figure 1 ). We found that our method required substantially fewer training samples (∼1e4 samples compared to PLDSID's ∼1e5) to reach ideal (i.e., ground truth) prediction (figure 1a ). Similarly, our method more accurately identified the shared dynamical modes compared to PLDSID even when methods had increasingly more training samples (figure 1b ). In our second simulation experiment, we studied the effect of latent state dimension on learning. We generated 50 systems with fixed dimensions for shared and total latent states given by n 1 = 4 and n x = 16, respectively. We swept the learned latent state dimension from 1 to the true dimensionality of n x = 16, with the dimensionality of shared dynamics set to min(current n x , n 1 ). We found that our method identified the correct shared modes with very small errors using only 4 latent state dimensions; in contrast, PLDSID did not reach such low error rates even when using higher latent state dimensions (figure 1d ). In terms of predictive power, our method achieved close-to peak performance even when using as few as 4 latent states whereas PLDSID required much larger latent states dimensions, of around 16, to do so (figure 1c ). Taken together, these results show the power of PG-LDS-ID for performing dimensionality reduction on Poisson observations while prioritizing identification of shared dynamics with a secondary Gaussian data stream.As a demonstration on real data, we used our algorithm to model the shared dynamics between discrete population neural spiking activity and continuous arm movements in a publicly available NHP dataset from the Sabes lab (O'Doherty et al., 2017) . The dataset is of a NHP moving a 2Dcursor in a virtual reality environment based on fingertip position. We use the 2D cursor position and velocity as the continuous observations z. We removed channels that had average firing rates less than 0.5 Hz or greater than 100 Hz. Similar to Lawlor et al. ( 2018), we also removed channels that were correlated with other channels using a correlation coefficient threshold of 0.4. For all methods we used 50ms binned multi-unit spike counts for the discrete observations y. We evaluated decoding performance of learned models using five-fold cross validation across six recording sessions. We performed cross-validation using randomly-selected, non-overlapping subsets of 15 channels (n y = 15) 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 i y , continuous horizon i z , and time lag, which specifies how much the neural time-series should be lagged to time-align with the corresponding behavioral time-series (Moran & Schwartz, 1999; Shoham et al., 2005; Pandarinath et al., 2018) . We swept i y values of 5 and 10 time bins, i z values of 10, 20, 22, 25, 28, and 30 time bins; and lag values of 0, 2, 5, 8, and 10 time bins.To train PG-LDS-ID, we use the shared dynamics dimensionality of n 1 = min(current n x , 8). We chose a maximum n 1 of 8 because behavior decoding roughly plateaued at this dimension.Compared with PLDSID, our method learned models that led to better behavioral decoding at all latent state dimensions (figure 2a ) and achieved a higher behavior decoding at the maximum latent state dimension. 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. residual latent states in Poisson observations. Interestingly, despite the focus on learning the shared latent states in the first stage, PG-LDS-ID was also able to extract the residual latent states in Poisson observations because of its second stage. This led to PG-LDS-SID performing similarly to PLDSID in terms of peak neural self-prediction AUC while outperforming PLDSID in terms of peak behavior decoding (figure 2c ). Indeed, even with the inclusion of just two additional latent states to model residual Poisson dynamics (n 2 = 2, n x = 10), neural self-prediction was comparable to models learned by PLDSID (figure 2b ). Taken together, our method was extensible to real data and helped boost decoding performance, especially in low-dimensional latent regimes, by better identifying shared dynamics between Poisson and Gaussian observations. In appendix A.8 we also include preliminary results comparing against PLDS models fit using EM on a subset of this dataset.PG-LDS-ID, similar to other SSID methods, uses a time-invariant model which may not be suitable if the data exhibits non-stationarity, e.g., in chronic neural recordings. In such cases one would need to intermittently refit the model or develop adaptive extensions (Ahmadipour et al., 2021) . Moreover, as with other covariance-based SSID methods, our method may be sensitive to 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 figure 1a-b . Further, it may be possible that SSID methods fail to learn a valid set of parameters corresponding to a positive-definite covariance sequence (Van Overschee & De Moor, 1996; Katayama, 2005) or they may learn unstable state dynamics, meaning A has some eigenvalues with magnitude greater than 1 (see appendix A.6). These issues can arise due to errors in the empirical estimates of the covariances and because these methods do not explicitly impose stability constraints on model parameters. Future work may consider incorporating techniques from control theory, such as mode stabilization and covariance matching, to help mitigate these limitations (Maciejowski, 1995; Lindquist & Picci, 1996; Byrnes et al., 1998; Alkire & Vandenberghe, 2002) . Finally, our modeling approach can only provide an approximation of nonlinear dynamics/systems within the class of generalized-linear models, which have been shown to well-approximate nonlinear data in many applications, including modeling of neural and behavioral data. 
discussion: 
 We developed a novel analytical two-staged subspace identification algorithm termed PG-LDS-ID for modeling Poisson data streams while dissociating the dynamics shared with Gaussian data streams. Using simulations and real NHP data, we demonstrated that our method successfully achieves this new capability and thus, compared to existing Poisson SSID methods, more accurately identifies Poisson dynamics that are shared with Gaussian observations. Furthermore, this capability allows our method to improve decoding performance despite using lower-dimensional latent states and requiring a fewer number of training samples. Although we specifically focused on modeling Gaussian and Poisson observations, our algorithm can be extended to alternate distributions described with generalized-linear models. Our algorithm only requires the second-order moments after moment conversion (see equations ( 2), (3), ( 7), ( 8)). Because the moment conversion algorithm can be modified for the desired link function in generalized-linear models, as explained by Buesing et al. (2012) , we can combine our method with the appropriate moment conversion to extend it to other non-Gaussian and non-Poisson observation distributions. Due to the high-prevalence of generalized-linear models across various application domains (e.g., biomedical engineering, neuroscience, finance, etc.), our method can be a general tool for modeling shared and residual dynamics of joint data streams with distinct observation distributions. 
reproducibility statement: 
 We have taken a few steps to ensure reproducibility of the results reported here. First, we are sharing the implementation of our algorithm, as supplementary material, along with example simulated data and a tutorial IPython notebook to demonstrate usage. Second, we used a publicly available dataset (O'Doherty et al., 2017) that can be easily accessed by anyone interested in reproducing the results reported in section 4.2. Finally, to further aid in reproducing results, we have also outlined the preprocessing and analyses steps we have taken in section 4.2 and appendix A.5.2.A.1 DERIVATION Here we provide the derivation for a prioritized covariance-based subspace identification algorithm that learns a dynamical model of a predictor data stream while dissociating shared vs. residual latents with a secondary data stream. We define the following equivalent formulation for our dynamical model (equation ( 5)), where the block structure delineates shared and unshared latent statesEQUATIONwith parameters and noise terms defined as in sections 2.1 and 3.1.Before we present the derivation for PG-LDS-ID, we review a few steps in standard covariancebased SSID (section 2.2) that will help us in the derivation. First, it can be shown that the τ -th lag cross-covariance terms for r can be written in terms of model parameters as Λ rτ = Cov (r k+τ , r k ) = C r A τ -1 G, where G := Cov(x k+1 , r k ). Using this relationship, the Hankel matrix, H, can be expanded as (Van Overschee & De Moor, 1996; Katayama, 2005 )EQUATIONSecond, using a singular-value decomposition, the above Hankel matrix H can be decomposed into observability, Γ r , and controllability, ∆, matrices from which model parameters can be extracted (Van Overschee & De Moor, 1996; Katayama, 2005 )H SVD = Γ r ∆ =     C r C r A . . . C r A i-1     A i-1 G • • • AG G . (16) A.1.2 PRIORITIZED COVARIANCE-BASED SUBSPACE IDENTIFICATION: STAGE 1 DERIVATIONIn 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 x(1)EQUATIONr , C z , b). First, it can be shown that the τ -th lag cross-covariance between z and r can be written in terms of model parameters asΛ zrτ = Cov (z k+τ , r k ) = C z A τ -1 G,where G is defined as before. Due to the block structure of equation ( 14), we have shown that G can be partitioned as G = G (1) G (2) (see section 3.1). As a result, ∆ (equation ( 16)) also has a block-partition formatEQUATION)and the cross-covariance term can be simplified asΛ zrτ = C (1) z 0 A 11 0 A 21 A 22 τ -1 G (1) G (2) = C (1) z A τ -1 11 G (1) .Henceforth we drop the superscript (1) on C(1)z , without loss of generality. The Hankel matrix between future continuous observations and past log-rates of the discrete observation can then be expanded asH zr = Cov(z f , r p ) =    Λ zri Λ zri-1 • • • Λ zr1 . . . . . . • • • . . . Λ zr2i-1 Λ zr2i-2 • • • Λ zri    =    C z A i-1 11 G (1) C z A i-2 11 G (1) • • • C z G (1) . . . . . . • • • . . . C z A 2i-2 11 G (1) C z A 2i-3 G (1) • • • C z A i-1 11 G (1) .    .A singular-value decomposition of H zr yields the observability matrix for z (i.e., Γ z ) and the controllability matrix ∆ (1) associated with the shared dynamicsEQUATION)At this point, C z can be read off the first n z rows of Γ z . The shared latent dynamics matrix A 11 can be learned by solving a least-squares problem based on the controllability matrix ∆ (1) (as introduced in section 3.2.1) 1) where ( 19)EQUATION∆ (1) := A i-1 11 G (1) • • • A 11 G (1) , ∆ (1) := A i-2 11 G (1) • • • G (1) ,which has the following closed-form solution:A 11 = ∆ (1) (∆ (1) ) † .To extract C r , we first note that the Hankel expansion in equation ( 16) can be simplified due to the block-structure of ∆ and Γ r , which is defined asΓ r =     C r C r A . . . C r A i-1     =       C (1) r C (2) r C (1)r A C(2) r A . . . . . .(1)EQUATION)using the fact that C r = C(1) r(2) r. Then, we explicitly separate H = Γ r ∆ (equation ( 16)) into two parts based on its singular-value decomposition asH = UΣV T = Γr (UΣ 1/2 ) ∆ (Σ 1/2 V T ) (a) = U (1) U (2) Σ 1/2 (1) 0 0 Σ 1/2 (2) Σ 1/2 (1) 0 0 Σ 1/2 (2)V T(1)V T(2)EQUATIONwhere simplification (a) is due to the block-partition structure of Γ r and ∆. Thus H can be written as the sum of "shared" and "residual" components (equation ( 21)). We can compute Γ(1) r asEQUATION)where we have used the orthonormal property of right singular vectors V to conclude ∆ (2) ∆ (1) † = 0. At this point, we can extract C(1) r by reading the top n y rows of Γ(1) r . Finally, b is learned directly during the moment transformation (section 2.1). This concludes the learning of all parameters associated with the shared dynamical subspace, i.e., (A 11 , C(1)r , C z , b). A.1.3 PRIORITIZED COVARIANCE-BASED SUBSPACE IDENTIFICATION: STAGE 2 DERIVATIONIn the second stage of our algorithm, our goal is to learn model parameters that describe the residual dynamics of r via the latent state x(2)EQUATIONr ). To learn these parameters, we first extract the residual component in equation ( 21), termed H (2) , by subtracting Γ(1) r ∆ (1) from H, and decompose it via a singular-value decomposition to get Γ(2) r and ∆ (2) asEQUATIONAt this point, we take C(2) r as the top n y rows of Γ(2) r and concatentate with C(1) r to complete C r . To complete the state dynamics matrix A, we refer back to the block-structure representation of the controllability matrix in equation ( 17)∆ (1) ∆ (2) = A 11 0 A 21 A 22 A i-2 G • • • A 11 0 A 21 A 22 G (1) G (2) G (1) G (2)from which we construct the following relationship∆ (1) ∆ (2) = A 11 0 A 21 A 22 ∆ (1) ∆ (2) (24)where ∆ and ∆ are defined as in equation ( 19). We can further isolate the residual state transitions as the solution to the following equation (taken from the second row of equation ( 24))EQUATIONwhich has the following closed-form least-squares solution:[A 21 A 22 ] = ∆ (2) ∆ † . The full state dynamics is the concatenation A = A 11 0 A 21 A 22. This concludes the learning of all parameters for the residual dynamics, i.e., ([EQUATIONr ).In sections A.1.1-A.1.3 we demonstrated how all model parameters can be extracted in two stages with prioritization, starting from the second-moments of the latent log-rates, r, and the Gaussian observations, z. Here we explain how we estimate these moments from the computable moments of y (the Poisson observations) and z, using equation ( 7) asΛ z fm rp n = Cov(z fm , y pn ) / µ yp n .Here we provide a sketch of the proof. Without loss of generality, assume z and r are stationary with a mean of 0 (e.g., demeaned during preprocessing). We can compute the covariance of any two elements j and k of vectors z f and y p respectively asCov(z fj , y p k ) = E z fj y p k = E E z fj y p k |r p k (a) = E E z fj |r p k E [y f k |r p k ] = E E z fj |r p k exp(r p k )where (a) is because z fj and y p k are independent when conditioned on latent log-rate r p k . Next, we use the fact that z f and r p are jointly Gaussian random processes and, as a result, the mean of the conditional distribution,E[z fj |r p k ], is equal to Λ z f j rp k Λ -1 rp kk r p k (i.e., the linear least-square estimate of z fj using r p k ). The last step is to compute the expectationE Λ z f j rp k Λ -1 rp kk r p k exp(r p k ) = Λ z f j rp k µ yp k = Cov(z fj , y p k )which, after rearranging terms, yieldsΛ z fm rp n = Cov(z fm , y pn ) / µ yp n .We note that the final equation is equivalent to a derivation provided by Buesing et al. (2012) as their supplementary equation ( 6) to compute cross-covariances between Poisson observations and Gaussian inputs, 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 (3) in section 2.1.For ease of exposition, the derivation in section A.1.2 was provided for a cross-term Hankel matrix H zr that was formed with equal horizons for z and r asH zr := Cov(z f , r p ) =     Λ zri Λ zri-1 • • • Λ zr1 Λ zri+1 Λ zri • • • Λ zr2 . . . . . . • • • . . . Λ zr2i-1 Λ zr2i-2 • • • Λ zri     , z f :=    z i . . . z 2i-1    , r p :=    r 0 . . . r i-1    .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 (Van Overschee & De Moor, 1996; Katayama, 2005) , i.e., n 1 = rank(H zr ) per equation ( 18) and n x = rank(H) per equation ( 16). 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 1 singular values in H zr ) 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 ( 16), (18), and (23), 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(H zr ) = min(i × n y , i × n z ) and rank(H) = i × n y . 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 PG-LDS-ID, we generalize the Hankel matrix H zr to support different horizon values for each of the observations, i z and i y , such thatH zr =      Λ zri z Λ zri z -1 • • • Λ zri z -iy +1 Λ zri z +1 Λ zri z • • • Λ zri z -iy +2 . . . . . . • • • . . . Λ zr2i z -1 Λ zr2i z -2 • • • Λ zr2i z -iy      with z f :=    z iz . . . z 2iz-1    , r p :=    r 0 . . . r iy-1    .The discrete observation horizon i y is also used when forming the Hankel matrix H, per equation ( 15). The additional flexibility gained from having different horizon values is especially critical in scenarios wherein the dimensionalities of z and y are very different, such as in the case of our NHP analysis where n z = 4 and n y = 15. We select the final horizons i z and i y via an inner cross-validation based on which values achieve the best accuracy in the training data.Here we explain how the blocked formulation in equation ( 6) can be assumed without loss of generality. The latent states in our model describe the primary data stream (y k , e.g., Poisson spiking activity), with a subset also explaining the secondary stream (z k , e.g., behavior). Formally, we define the true dimensionality of the shared states (denoted by n 1 ) based on the rank of the observability matrix for the pair (A, C z ). 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 n 1 dimensional latent subspace that is observable via z k as the first few dimensions of the latent space, thus giving the block-structured formulation of equation ( 6). This can be seen by applying Theorem 3.8 from Katayama (2005) to the first two lines of equation ( 5). Thus, the blocked formulation of equation ( 6) is equivalent to the formulation from ( 5) and we can aim to learn our model in the form of equation ( 6) without any loss of generality. Moreover, note that this blocked formulation also covers the special case of a non-blocked formulation when n 1 = n x , that is when all latent states contribute to both data streams and the observability matrix is full-rank. In this case, the top-left-block of A grows to cover the whole A, and thus no zero-filled upper-right block would remain. The algorithm would still work in this special case by simply only applying the first stage of learning. However, within the application of modeling neural and behavioral data, we typically expect a minority of the neural dynamics to be related to a particular behavior of interest and so we expect the most appropriate n 1 to be smaller than n x .Hyperparameters n 1 and n x denote the number of shared vs. total latent state dimensions (equation ( 6)). When modeling real data, one can estimate the most appropriate values for these hyperparameters for the data using the following procedure:1. Sweep over values of n 1 , increasing n 1 while keeping n x = n 1 (i.e., using stage 1 only to learn). Quantify the prediction of the secondary data stream (e.g., behavior) in each case to find the n 1 at which the prediction plateaus or reaches a peak. This value gives the appropriate n 1 . Alternatively, n 1 can be estimated as the number of non-zero (or nonnegligible) singular values of the Hankel matrix H zr (equations ( 8) and ( 9)). 2. Using the selected n 1 from above, sweep over values of n x , starting from n 1 and increasing the latent state dimension. Quantify the self-prediction of the predictor data stream (e.g., spiking activity) in each case, and find the n x at which the self-prediction reaches a peak.The model parameters stated in section 3.1 correspond to the parameters required by the pointprocess filter (Eden et al., 2004) for state estimation (see appendix A.4). However, if desired, under Gaussian assumptions the noise covariance term for the secondary data stream, Λ ϵ , can be learned by computing the covariance of the prediction residuals as E[(ẑ k -z k )(ẑ k -z k ) T ],where ẑk denotes the predicted value of z and the expectation denotes the empirical average across all time samples. When ϵ k is not white, the behavior prediction residuals under Gaussian assumptions can be computed in the same way (i.e., ẑkz k ) and modeled using Gaussian SSID.The Gaussian process' noise has an effect on learning similar to that of the residual dynamics present in the primary (i.e., Poisson) data stream: it reduces the overall signal-to-noise ratio of the shared dynamics present in both observation data streams. However, with increasingly more training samples and improved estimates of the second-order moments (i.e., covariances), the algorithm can become more robust to the impact of Gaussian observation noise.For model evaluations we chose to predict the continuous Gaussian observations from the discrete Poisson observations, which is a common use-case in neuroscience (Koyama et al., 2010; Macke et al., 2015; Lu et al., 2021) . Once model parameters are learned in a training set, either with PLDSID or PG-LDS-ID, we can use the learned parameters to construct the Poisson point-process filter (PPF) (Eden et al., 2004) and estimate the latent states in a test set. Note, using the PPF for state estimation is only possible if noise statistics are valid (section 3.2.3). We denote the one-step ahead latent state prediction of x k using all samples of y k up to time k -1 by xk|k-1 . These state estimates can be used to predict the continuous observations as C z xk|k-1 . To learn a C z parameter for PLDSID, we first estimate the latent states in the training data using a PPF and then fit a linear regression (scikit-learn) from the latent states to z k , i.e., C z = Z XT ( X XT ) † , where columns of Z and X contain z k and xk|k-1 for all training timepoints k (Pedregosa et al., 2011) . To make the methods more comparable, we use the same approach to refit the C z learned by PG-LDS-ID. We quantify the decoding performance using correlation coefficient (CC). We also assessed the one-step ahead self-prediction of Poisson observations using the predicted latent states. This was quantified with the area under the curve (AUC) of the receiver operating characteristic (appendix A.5.3).A.5 EXPERIMENTAL DETAILS A.5.1 SIMULATIONS For our synthetic data in section 4.1, we simulated Poisson-Gaussian observations from random models as per equation ( 5). We randomly selected the latent state dimension n x , the shared dimension n 1 , and the observation dimensions n y and n z with uniform probability from the following ranges: 1 ≤ n x ≤ 10, 20 ≤ n y ≤ 30, 5 ≤ n z ≤ 10, and 1 ≤ n 1 ≤ n x . Using these dimensions we generated random model parameters Θ = (A, C r , C z , b, Q). We constrained the complex eigenvalues (i.e., 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 (Churchland et al., 2012; Song et al., 2022) . All simulations were on a 10 ms timescale with a baseline log rate, b, randomly selected within [0.5, 15] Hz. Observation matrix C r was scaled to achieve a desired per-dimension maximum firing rate such that max 1≤k≤N exp(r k ) ∈ [25, 65] Hz. This was to ensure a realistic range of firing rate and sufficient modulation depth for all dimensions of the simulated Poisson point-process. Q was randomly generated to be a positive definite matrix and C z was generated to hit a target signal-to-noise, defined as the variance associated with latent states normalized by observation noise variance (C z Λ x C T z )/(Λ ϵ ). Target SNR values were randomly generated as 10 α with α uniformly distributed between [0, 2]. For every simulated model, the colored noise for the Gaussian process, ϵ k , was taken as the output of a 4-dimensional latent linear dynamical system with random parameters that were generated similarly. By using a general colored noise, we can simulate dynamics present in the continuous modality that are unshared with the discrete modality. To evaluate how well our algorithm modeled neural dynamics, we computed the one-step ahead selfprediction performance as quantified by AUC (section A.4). Our goal was to validate if our model could, when using all past neural observations, accurately predict the occurrence of spikes versus no spikes in a given time step. Since all self-predictions are made using the recursive point-process filter estimates of the latent states (see appendix A.4), we computed the probability of a spiking event for the m-th dimension of y at time k, conditioned on all observations y 1:k-1 , asP (y m k > 0|y 1:k-1 ) = x k p(y m k > 0 | y 1:k-1 , x k )p(x k | y 1:k-1 ) (a) = E x k | y 1:k-1 [p(y m k > 0 | x k )] (b) = E x k | y 1:k-1 [1 -exp(exp(r m k )) | x k ] (c) ≈ E x k | y 1:k-1 [exp(r m k ) | x k ] (d) = exp rm k|k-1 + 1 2 Λ rmmwhere in (a) we simplify using y k 's conditional independence from the past y 1:k-1 , in (b) we simplify based on y k | x k ∼ Poisson(exp(r k )), in (c) we use the Taylor series approximation of exp(exp(r k )) for small exp(r k ), and (d) is simply the mean of a log-normal random variable. Note that rk = C r xk|k-1 + b and Λ r = C r Λ xk|k-1 C T r , where xk|k-1 is the current estimate for the state and Λ xk|k-1 the state-prediction covariance (appendix A.4).Subspace identification methods generally only converge to the correct system parameters asymptotically (see figure 1a-b ), as the empirically estimated covariances also converge to their true values (Van Overschee & De Moor, 1996) . 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. In simulations (figure 1 ), we excluded learned models that were unstable from the reported mean performances, reflecting in the reduced number of samples in the standard error of the mean (s.e.m). For training set sizes typical of neuroscience datasets, the occurrences of unstable models was rare, with only 2 unstable systems for 1e5 training samples and no unstable systems for 1e6 training samples.We measured the learning time of our method and, additionally, the inference time associated with using a point-process filter for state estimation. Using one session of NHP data (section 4.2), we repeatedly trained on 25 distinct time-series datasets. Each dataset consisted of a 6097-by-15 matrix (timesteps-by-features) of Poisson observations and a 6097-by-4 matrix of Gaussian observations. The training time of our algorithm averaged across 25 trials was 0.33s (including 0.072s spent on the convex optimization problem outlined in 3.2.3). Further, we also measured the inference time for our 1524-sample testset to be 0.33s (i.e., approximately 0.2ms per 50ms timestep).Most of the computational cost of our algorithm is involved in the matrix operations associated with 1) computing the necessary covariance/Hankel 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.1) and a future-past Gaussian-Poisson Hankel matrix (section 3.2.1). Both of these empirical estimates 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 here, wherein n y = 15 and i y = 10 (horizon). The computed square Poisson-Poisson covariance matrix was of dimension 2 * n y * i y = 2 * 10 * 15 and was the result of a matrix multiplication between two matrices of dimension (2*10*15)-by-6078, where 6078 = timesteps -2 * i y +1. 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. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: * The multi-step nature of the algorithm may introduce significant computational overhead, potentially limiting its scalability.
* The method assumes time-invariant dynamics, which may not hold for all neural data, potentially affecting its generalizability.
* 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.
questions: 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.
* 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: limitations: To my knowledge and understanding, there is no potential negative societal impact of this work. Other limitations please see “Weaknesses”.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
questions: questions: See Weaknesses.
limitations: limitations: See Weaknesses.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

I 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.
The 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.

The 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?

There 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.
questions: 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? 

I 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? 

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

Did the matrices $Q$ described in line 998 conform with assumption A.1 ? 

In A.7.1, why describe the norm as Froebenius rather than Euclidean when it's about vectors? 

Choosing red and green as important colors for the figures is not ideal for colorblind people.
limitations: limitations: The authors discuss limitations of their work in a dedicated section.

The 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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:
$ 
TA^{-1}T = \begin{bmatrix} A_{11} & A_{12} \\\\ 0 & A_{22} \end{bmatrix}, CT = \begin{bmatrix}0& C_2\end{bmatrix}
$, 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.  
- 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.

Minor

- 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.
questions: 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?   
- Why does PGLDM always perform worse at self-prediction compared to other baselines (Figure 1 (b), 2(b), 3(b))?
- 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: 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.

Comments: comment: I thank the authors for their response and the clarifications. I will keep my rating as is.
We thank the reviewer for following up. We are glad the assumptions are clearer now and that the reviewer finds them interesting. We will more explicitly discuss them in the manuscript, as suggested. The reviewer is correct that the reasoning behind the assumption $A_{31}=0$ is primarily due to designation of one time-series as primary and the other as secondary. Moreover, the reviewer raises an interesting point about the bio-physical interpretation of these assumptions. We agree with their interpretation that the primary/secondary designation implies, roughly, a direction of causality. The reviewer specifically highlights two scenarios: unidirectional interaction – with either known or unknown direction – and bidirectional interaction (e.g., feedforward and feedback). The designation of primary/secondary time-series has a bio-physical basis in applications wherein there is clear unidirectional interaction, consistent with the reviewer’s intuition. For example, one important use-case of our method is dissociating behaviorally-relevant neural dynamics. In this case, it can make sense to designate the neural activity as the primary time-series and behavior as the secondary time-series, with the reasonable bio-physical assumption that the brain drives behavior, in a unidirectional manner. Regarding the distinct use-case of modeling shared dynamics between brain regions, we agree with the reviewer that this designation makes most bio-physical sense in applications where the direction is from one brain region to the other, rather than being simultaneously in both directions (both feedforward and feedback). If there is a unidirectional interaction and the direction is known, the primary will become the upstream region activity and the secondary the downstream region activity. If the direction is unknown, one may be able to build models for each alternative designation and compare them according to the desired goodness-of-fit criterion. Finally, although the reviewer is correct that our current model formulation would not be able to handle bidirectional communication simultaneously, it can still be used to model each direction separately (i.e., model both possible designations of primary/secondary time-series separately). Regardless, we agree with the reviewer that being able to model both communication directions simultaneously could be even more general and will be an interesting direction for future investigation.

Overall, we thank the reviewer for their very helpful feedback; we enjoyed these thoughtful discussions. We will definitely include explicit discussion regarding our model assumptions in the manuscript as well as their potential bio-physical implications – as the reviewer suggested.
I now better understand the authors' assumptions. I think the work will benefit by more explicitly discussing these assumptions. Regardless, these explanations do not alleviate my concerns over the bio-physical interpretation of these assumptions. One of the attractive property of linear dynamical models over non-linear models is their interpretability. The reasoning behind the assumption $A_{31}=0$ does not seem to stem from a bio-physical origin but as I understand (please clarify if I am wrong) from the fact that their is a primary time-series (which is designated by the "experimentalist/data scientist") and an additional "secondary" time-series which should not be able to predict the private dynamics ($x^{(2)}$) of the "primary" time-series. The designation of the "primary" and "secondary" time-series by itself sort of assumes a direction of causality/communication (I am using the words a bit vaguely) which in many cases might not be known a priori (many dynamical systems are indeed used for discovering this direction, see numerous works from the field of dynamic causal modeling). While this designation of time-series (e.g., from different brain regions) being primary and secondary might be perfectly reasonable for some applications (e.g., where there is a brain region sending information to another region without the other brain region providing feedback), there are also many cases where such an assumption might not be reasonable (e.g., in cases where brain areas communicate in both directions). This makes the model assumptions more restrictive. I still think the model assumption explored by the authors are interesting but can be restrictive.
We thank the reviewer for reading the rebuttal and for their thoughtful follow-up. We would like to clarify a few points:
1. We acknowledge that assuming $A_{31}=0$ simplifies the derivation because it allows the optional third stage (i.e., extraction of the secondary time-series’ private dynamics) to operate independently of the first stage (i.e., extraction of the shared dynamics). We initially realized this relationship between $A_{31}$ and stage 3 while working on the derivation, and determined that this scenario (i.e., $A_{31} \neq 0$) would not be of interest in our use-case since having a non-zero $A_{31}$ also meant that y (or r) would be predictive of private t (or z) dynamics (i.e., $x^{(3)}$). We will, as per the reviewer’s point, more explicitly state in the manuscript that the $A_{31}=0$ model assumption simplifies the algorithm derivation.

2. We do not claim that $A_{31}$ and $A_{21}$ are fundamentally different in the sense of information leakage, but we clarify that their roles in our problem formulation are asymmetric for the following reason. In our formulation, y/r is designated as the primary time-series (predictor) and t/z as the secondary time-series (predicted). Given this designation, we always exclusively use y/r as the predictor and t/z as the predicted; this means that a non-zero $A_{21}$ does not allow t/z to directly predict future private dynamics of y/r (i.e., $x^{(2)}$). However, having a non-zero $A_{31}$ allows y/r to directly predict future private dynamics of t/z. This distinction between $A_{21}$ and $A_{31}$ arises because during learning we only ever compute the lagged cross-correlations $\Lambda_{zr_{\tau}} = \mathrm{Cov}(z_{k+\tau}, r_{k})$ (eq 8) and not $\Lambda_{rz_{\tau}} = \mathrm{Cov}(r_{k+\tau}, z_{k})$ -- as per the primary/secondary designation. Regardless, our solution can handle the symmetric case when both $A_{31}$ and $A_{21}$ are zero as well, see next item.

3. Regarding the $A_{21} = 0$ case, we agree with the reviewer that this is a very valid formulation as well and can be of interest to users. Our solution can learn a model with this block-diagonal $A$ structure with both $A_{21} = 0$ and $A_{31} = 0$ with a minor modification to the second stage, as we noted in our rebuttal. We have now also tested the block-diagonal $A$ case on a single session of the primate dataset in Fig 2 and confirmed that the block-diagonal $A$ formulation also identified the shared dynamics more accurately than baselines and performed comparably to the original non-block-diagonal $A$ formulation (eq 7). We will include discussion on the merits of this block-diagonal $A$ formulation in our manuscript, based on the reviewer’s excellent suggestion, so that users can choose whichever option they prefer in their application (a non-zero $A_{21}$ or $A_{21}=0$).
I have read the rebuttal. I am still confused about the reasoning provided by the authors regarding the structure of A matrix. I am not entirely sure about the notion of private and shared dynamics authors are assuming. While I am okay with the authors' reasoning that $A_{12}=A_{13}=A_{23}=A_{32}=0$. The reasoning for $A_{31}=0$ seems very forced. I do not understand authors' reasoning that since they only use $y_k$ during inference, it is crucial that $y_k$ should not leak information about the other time-series $z_k$ through the shared dynamics and $A_{31}$ whereas since $t_k$ is not used during inference, it is okay for $t_k$ to leak information about the private dynamics through the shared dynamics. It seems very unlikely that neuronal dynamics in the brain would be affected by which time-series is being observed in the experiment. Either go all the way in assuming that shared dynamics and private dynamics do not leak information about each other or use the general structure. In my opinion, it seems more likely that $A_{31}=0$ is a simplifying assumption which is a *fine* assumption to make, as empirical evidence shows that the method work. I am keeping my original score.