Spectral Learning of Shared Dynamics Between Generalized-Linear Processes

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.

Generality claim on block form structure described in eq (7)

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

1. 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$.
2. 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$.
3. 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.

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

Structure of the proposed $A$ matrix

We thank the reviewer for their great question. We clarify a few points:

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

PGLDM self-prediction performance

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

V2 decoding from V1; V1 self-prediction

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

Additional comments

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

References:

[1] Tohru Katayama. Subspace Methods for System Identification.

[2] Semedo, J.D., Jasper, A.I., Zandvakili, A. et al. Feedforward and feedback interactions between visual cortical areas use different population activity patterns.

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.

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.

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.

Feedback on Def 3.1, section A.2, & assumption A.1

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

accumulation of errors in the presence of estimation noise

We thank the reviewer for their great question. There are two aspects of the learning we wish to comment on.

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

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

clarification on the first experiment

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

bad performance of EM in Fig 1

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

Structure of synthetic $Q$ matrices

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

Additional comments

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

References:

[1] Linderman lab SSM library

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.

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

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

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

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.

The method assumes time-invariant dynamics, which may not hold for all neural data, potentially affecting its generalizability.

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

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.

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

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

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

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

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

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

References:

[1] Parima Ahmadipour, Yuxiao Yang, Edward F. Chang, and Maryam M. Shanechi. Adaptive tracking of human ECoG network dynamics.