Identifying Feedforward and Feedback Controllable Subspaces of Neural Population Dynamics

## Abstract

Abstract There is overwhelming evidence that cognition, perception, and action rely on feedback control. However, if and how neural population dynamics are amenable to different control strategies is poorly understood, in large part because machine learning methods to directly assess controllability in neural population dynamics are lacking. To address this gap, we developed a novel dimensionality reduction method, Feedback Controllability Components Analysis (FCCA), that identifies subspaces of linear dynamical systems that are most feedback controllable based on a new measure of feedback controllability. We further show that PCA identifies subspaces of linear dynamical systems that maximize a measure of feedforward controllability. As such, FCCA and PCA are data-driven methods to identify subspaces of neural population data (approximated as linear dynamical systems) that are most feedback and feedforward controllable respectively, and are thus natural contrasts for hypothesis testing. We developed new theory that proves that non-normality of underlying dynamics determines the divergence between FCCA and PCA solutions, and confirmed this in numerical simulations. Applying FCCA to diverse neural population recordings, we find that feedback controllable dynamics are geometrically distinct from PCA subspaces and are better predictors of animal behavior. Our methods provide a novel approach towards analyzing neural population dynamics from a control theoretic perspective, and indicate that feedback controllable subspaces are important for behavior.

## 1 Introduction

Feedback control has long been recognized to be central to brain function [Wiener, , Conant and Ashby, ]. Prior work has established that, at the behavioral level, motor coordination [Todorov and Jordan, ], speech production [Houde and Nagarajan, ], perception [Rao and Ballard, ], and navigation [Pezzulo and Cisek, , Friston et al., ] can be accounted for by models of optimal feedback control. Advances in the ability to simultaneously record from large number of neurons have further revealed that the brain computes through population dynamics [Vyas et al., ]. Nonetheless, whether neural population dynamics have particular subspaces that are more or less amenable to feedback control is unknown. Addressing this gap requires the development of novel methods to identify components of population activity relevant for control.

The cost incurred in controlling a dynamical system is referred to as its controllability. Existing measures of controllability center around the energy (in terms of the norm of the control signal) that must be expended to steer the system state. These measures are calculated from the controllability Gramian of the (linearized) system dynamics. Controllability is an intrinsic feature of the dynamical system itself, and may be estimated from measurements of system dynamics without reference to the specific inputs to the system Pasqualetti et al. []. Network controllability analyses have delivered insights into the organization of proteomic networks [Vinayagam et al., ], human functional and structural brain networks [Medaglia et al., , Tang and Bassett, , Kim et al., , Gu et al., ], and the connectome of *C Elegans* [Yan et al., ]. However, prior work in network controllability has exclusively focused open loop, or feedforward, controllability in the context of extracted networks, and not measures of closed loop, or feedback, controllability in the context of observed dynamics of data. Indeed, methods to asses feedback controllability from observations of the dynamics of neural populations are nascent.

Here, we developed dimensionality reduction methods that can be applied to neural population data that maximize the feedforward and feedback controllability of extracted subspaces. We first identify a correspondence between Principal Components Analysis (PCA) and the volume of state space reachable by feedforward control in linear dynamical systems [Pasqualetti et al., ]–this provides a control-theoretic interpretation to PCA extracted subspaces. We then present Feedback Controllable Components Analysis (FCCA), a linear dimensionality reduction method to identify feedback controllable subspaces of high dimensional dynamical systems based on a novel measure of feedback controllability. We show that the FCCA objective function can be applied to data using only its second order statistics, bypassing the need for prior system identification and making the method easily applicable to high dimensional neural population recordings.

Through theory and numerical simulations, we show that the degree of non-normality of the underlying dynamical system [Trefethen and Embree, ] determines the degree of divergence between PCA and FCCA solutions. In the brain, the postsynaptic effect of every neuron is constrained to be either excitatory or inhibitory by Dale’s Law. This structure implies that linearized dynamics within cortical circuits are necessarily non-normal [Murphy and Miller, ]. Prior work has highlighted the capacity of non-normal dynamical systems to retain memory of inputs [Ganguli et al., ] and transmit information [Baggio and Zampieri, ]. Our results show that non-normality also plays a fundamental role in shaping the controllability of neural systems. Finally, we applied FCCA to diverse neural recordings and demonstrate that feedback controllable subspaces are better predictors of behavior than PCA subspaces (despite both being linear), and that the two subspaces are geometrically distinct.

## 2 Controllable subspaces of linear dynamical systems

Here, we provide detailed derivations of our data-driven measures of controllability. We first discuss the natural cost function to measure feedforward controllabiity (eq. 4) and highlight its correspondence to PCA. Next, we present the analogous measure for feedback controllabiity (eq. 7), and how it may be estimated *implicitly* (i.e., without explicit model fitting) from the observed second order statistics of data (eq. 11). We provide rationale for this cost function as measuring the complexity of the feedback controller required to regulate the observed neural population dynamics.

We consider linear dynamical systems of the form:

$$ $\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+Bu(t)\quad y(t)=Cx(t)$ (1) $$

where $x(t)\in\mathbb{R}^{N}$ is the neural state (i.e., the vector of neuronal activity, not a latent variable) and $u(t)$ is an external control input. $A\in\mathbb{R}^{N\times N}$ is the dynamics matrix encoding the effective first order dynamics between neurons. $B\in\mathbb{R}^{N\times p}$ describes how inputs drive the neural state, and $C\in\mathbb{R}^{d\times N},d<<N$ is a readout matrix projecting the neural dynamics to a lower dimensional space. The input-output behavior (i.e., the mapping from $u(t)$ to $y(t)$) can equivalently be represented in the Laplace domain using the transfer function $G(s)=C(sI-A)^{-1}B$ Kailath [].

Consider an invertible linear transformation of the state variable $x\to Tx$. Under such a state-space transformation, the input-output behavior of the system 1 is left unchanged as the state space matrices transform as $(A,B,C)\to(TAT^{-1},TB,CT^{-1})$. This implies that there are many possible choices of $(A,B,C)$ matrices, referred to as realizations, that give rise to the same transfer function $G(s)$. A minimal realization contains the fewest number of state variables (i.e., $A$ has the smallest dimension) amongst all realizations. Measures of controllabity that are *intrinsic* to the dynamical system should be invariant across all realizations. We will show that our measures of feedforward and feedback controllabillity exhibit this property.

Throughout, we will assume that the observed data obeys following underlying state dynamics:

$$ $\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+Bdw(t);\quad dw(t)\sim\mathcal{N}(0,1);\quad y(t)=Cx(t)$ (2) $$

Compared to eq. 1, $u(t)$ has been replaced by temporally white noise $dw(t)$, a reasonable assumption given that input signals are unmeasured in neural recordings. Our metrics of controllability rely only on observing the linear dynamics under this latent, stochastic excitation.

### 2.1 Principal Components Analysis Eigenvalues Measure Feedforward Controllability

A categorical definition of controllability for a dynamical system is that for any desired trajectory from initial state to final state, there exists a control signal $u(t)$ that could be applied to the system to guide it through this trajectory. For a (stable) linear dynamical system, a necessary and sufficient condition for this to hold is that the controllability Gramian, $\Pi$, has full rank. $\Pi$ is obtained from the state space parameters through the solution of the Lyapunov equation:

$$ $\displaystyle A\Pi+\Pi A^{\top}=-BB^{\top}\quad\Pi=\int_{0}^{\infty}dt\;e^{At} BB^{\top}e^{A^{\top}t}$ (3) $$

The rank condition on $\Pi$ as a definition of controllability, while canonical [Kailath, ], is an all or-nothing designation; either all directions in state space can be reached by control signals, or they cannot. Furthermore, this definition does not take into account the energy required to achieve the desired transition. While certain directions in state space may in principle be reachable, the energy required to push the system in those directions may be prohibitive.

Thus, given that the system is controllable, we can ask a more refined question: what is the energetic effort required to control different directions of state space? The energy required for control is measured by the norm of the input signal $u(t)$. It can be shown [Pasqualetti et al., ] that to reach states that lie along the eigenvectors of $\Pi$, the minimal energy is proportional to the inverse of the corresponding eigenvalues of $\Pi$. Directions of state space that have large projections along eigenvectors of $\Pi$ with small eigenvalues are therefore harder to control. For a unit-norm input signal, the volume of reachable state space is proportional to the determinant of $\Pi$ [Summers et al., ].

The above intuition can be encoded into the objective function of a dimensionality reduction problem: for a fixed-norm input signal, find $C$ that maximizes the reachable volume within the subspace. This volume is measured by the determinant of $C\Pi C^{\top}$. Identifying subspaces of maximum feedforward controllability is then posed as the following optimization problem:

$$ $\displaystyle\text{argmax}_{C}\log\det C\Pi C^{\top}\;\;|\;\;C\in\mathbb{R}^{d \times N},CC^{\top}=I_{d}$ (4) $$

Observe that under state space transformations, $\Pi$ maps to $T\Pi T^{\top}$, whereas $C$ maps to $CT^{-1}$. Hence, as desired, eq. 4 is invariant to state space transformations and thus an intrinsic property of the dynamical system. We include the constraint $CC^{\top}=I_{d}$ to ensure the optimization problem is well-posed. Without it, one could, for example, multiply $C$ by a constant and increase the objective function. We can assess this objective function from data generated by eq. 2, as in this case the observed covariance of the data will coincide with the controllability Gramian [Mitra, , Kashima, ]. The solution of problem 4 coincides with that of PCA, as the optimal $C$ of fixed dimensionality $d$ has rows given by the top $d$ eigenvectors of $\Pi$ (see Theorem 2 on pg. 7 and Lemma 1 in the Appendix).

### 2.2 Linear Quadratic Gaussian Singular Values Measure Feedback Controllability

How does one quantify the feedback controllability of a system? The primary distinction between feedforward control and feedback control is that the latter utilizes observations of the state to synthesize subsequent control signals. Feedback control therefore involves two functional stages: filtering (i.e., estimation) of the underlying dynamical state ($x(t))$ from the available observations ($y(t)$) and construction of appropriate regulation (i.e., control) signals. For a linear dynamical system, state estimation is optimally accomplished by the Kalman filter, whereas state regulation is canonically achieved via linear quadratic regulation (LQR). It will be crucial in what follows to recall that the Kalman Filter is an efficient, recursive, Gaussian minimum mean square error (MMSE) estimate of $x(t)$ given observations $y(\tau)$ for $\tau\leq t$. These two functional stages optimally solve the following cost functions:

$$ $\displaystyle\text{Kalman Filter}:$ $\displaystyle\;\;\min_{p(x_{0}|y_{-T:0})}\lim_{T\to\infty}\text{Tr}\left( \mathbb{E}\left[(\mathbb{E}(x_{0}|y_{-T:0})-x_{0})(\mathbb{E}(x_{0}|y_{-T:0})- x_{0})^{\top}\right]\right)$ $\displaystyle\text{LQR}:$ $\displaystyle\;\;\min_{u\in L^{2}[0,\infty)}\;\;\lim_{T\to\infty}\mathbb{E} \left[\frac{1}{T}\int_{0}^{T}x^{\top}C^{\top}Cx+u^{\top}u\;dt\right]$ $$

where $y_{-T:0}$ denotes observations over the interval $[-T,0]$. The minima of these cost functions are obtained from the solutions of dual Riccati equations:

$$ $\displaystyle AQ+QA^{\top}+BB^{\top}-QC^{\top}CQ$ $\displaystyle=0$ (5) $\displaystyle A^{\top}P+PA+C^{\top}C-PBB^{\top}P$ $\displaystyle=0$ (6) $$

where

$$ $\displaystyle Q$ $\displaystyle=\min_{p(x_{0}|y_{-T:0})}\lim_{T\to\infty}\mathbb{E}\left[( \mathbb{E}(x_{0}|y_{-T:0})-x_{0})(\mathbb{E}(x_{0}|y_{-T:0})-x)^{\top}\right]$ $\displaystyle x_{0}^{\top}Px_{0}$ $\displaystyle=\min_{u\in L^{2}[0,\infty)}\left\{\lim_{T\to\infty}\mathbb{E} \left[\frac{1}{T}\int_{0}^{T}x^{\top}C^{\top}Cx+u^{\top}u\;dt\right],\;\;x(0)= x_{0}\right\}$ $$

Here, $Q$ is the covariance matrix of the estimation error, whereas $P$ encodes the regulation cost incurred for varying initial conditions ($x_{0}$). $\text{Tr}(P)$ is proportional to the average regulation cost over all unit norm initial conditions.

The solutions of the Riccati equations are not invariant under the invertible state transformation $x\mapsto Tx$. The filtering Riccati equation will transform as $Q\mapsto TQT^{\top}$ whereas $P$ will transform as $(T^{-1})^{\top}PT^{-1}$. As such, simply by defining new coordinates via $T$ we can shape the difficulty of filtering and regulating various directions of the state space. Therefore $Q$ and $P$ on their own are not suitable cost functions for measuring feedback controllability. However, the product $PQ$ undergoes a similarity transformation $PQ\to(T^{\top})^{-1}QPT^{\top}$. *Hence, the eigenvalues of $PQ$ are invariant under similarity transformations, and define an intrinsic measure of the feedback controllability of a system. Additionally, there exists a particular $T$ that diagonalizes $PQ$.* Following [Jonckheere and Silverman, ], we refer to the corresponding eigenvalues as the LQG (Linear Quadratic Gaussian) singular values. In this basis, the cost of filtering each direction of the state space equals the cost of regulating it. We formalize these statements by restating Theorem 1 from [Jonckheere and Silverman, ]:

###### Theorem 1 .

Let $(A,B,C)$ be a minimal realization of $G(s)$. Then, the eigenvalues of $QP$ are similarity invariant. Further, these eigenvalues are real and strictly positive. If $\mu_{1}^{2}\geq\mu_{2}^{2}\geq\mu_{N}^{2}>0$ denote the eigenvalues of $QP$ in decreasing order, then there exists
a state space transformation $T$, $(A,B,C)\to(TAT^{-1},TB,CT^{-1})\equiv(\tilde{A},\tilde{B},\tilde{C})$ such that:

$$ $Q=P=\text{diag}(\mu_{1},\mu_{2},...,\mu_{N})$ $$

The realization $(\tilde{A},\tilde{B},\tilde{C})$ will be called the closed-loop balanced realization.

###### Proof.

Let $Q=LL^{\top}$ be the Cholesky decomposition of $Q$ and let $L^{\top}PL$ have Singular Value Decomposition $U\Sigma^{2}U^{\top}$. Then, one can check $T=\Sigma^{1/2}U^{\top}L^{-1}$ provides the desired transformation.
∎

Hence, as an intrinsic measure of feedback controllability, we take the sum of the LQG singular values $\mu_{i}^{2}$, corresponding to the sum of the ensemble cost to filter and regulate each direction of the neural state space:

$$ $\displaystyle\text{Tr}(PQ)$ (7) $$

###### Theorem 1 .

###### Proof.

### 2.3 The Feedback Controllability Components Analysis Method.

We developed a novel dimensionality reduction method, Feedback Controllability Components Analysis (FCCA), that can be readily applied to observed data from typical systems neuroscience experiments. To do so, we construct estimators of the LQG singular values, and hence $\text{Tr}(PQ)$, directly from the autocorrelations of the observed neural firing rates. The FCCA objective function arises from the observation that causal and acausal Kalman filtering are also related via dual Riccati equations. We first show that through an appropriate variable transformation, we obtain a state variable $x_{b}(t)$ whose dynamics unfold backwards in time via the same dynamics matrix ($A$) which evolves $x(t)$ (the neural state) forwards in time. Once established, this enables us to use the error covariance matrix of Kalman filtering $x_{b}(t)$ as a stand-in for the cost of regulating $x(t)$.

In particular, given the state space realization of the forward time stochastic linear system in eq. 2
, the joint statistics of $(x(t),y(t))$ can equivalently be parameterized by a Markov model that evolves backwards in time [L. Ljung and T. Kailath, ]:

$$ $\displaystyle-\dot{x}_{b}(t)$ $\displaystyle=A_{b}x_{b}(t)+Bdw(t);\quad y=Cx_{b}(t)$ (8) $$

where $A_{b}=-A-BB^{\top}\Pi^{-1}=\Pi A^{\top}\Pi^{-1}$ and $\Pi=\mathbb{E}[x(t)x(t)^{\top}]$ is the solution of the Lyapunov equation (eq. 3)

Examination of eq. 5 and eq. 6 reveals that the filtering and LQR Riccati equations differ primarily in two respects. First, the dynamics matrix is transposed ($A\to A^{\top}$). Second, the inputs and outputs have been exchanged ($B\to C^{\top}$, $C\to B^{\top}$). To use the error covariance of state filtering as a stand-in for the state regulation cost, we therefore require that the corresponding acausal state dynamics (determined by $A_{b}$) respect these differences. To this end, consider the transformed state $x_{a}(t)=\Pi^{-1}x(t)$. Substituting $x(t)=\Pi x_{a}(t)$ and $A_{b}=\Pi A^{\top}\Pi^{-1}$ into the equations for the backward dynamics result in following dynamics for this *adjoint* state:

$$ $\displaystyle-\dot{x}_{a}(t)=A^{\top}x_{a}(t)+\Pi^{-1}Bdw(t)$ $$

Then, if we construct a readout of this transformed state $y_{a}(t)=C\Pi x_{a}(t)=Cx(t)$, the Riccati equation associated with Kalman filtering $x_{a}$, whose solution we denote $\tilde{P}$, takes on the form:

$$ $\displaystyle A^{\top}\tilde{P}+\tilde{P}A+\Pi^{-1}BB^{\top}\Pi^{-1}-\tilde{P} \Pi C^{\top}C\Pi\tilde{P}$ $\displaystyle=0$ (9) $\displaystyle A^{\top}P+PA+C^{\top}C-PBB^{\top}P$ $\displaystyle=0$ (eq 6 ) $$

We see that eq. 9 coincides with eq. 6 (reproduced for convenience) upon switching the inputs and outputs ($B\to C^{\top}$, $C\to B^{\top}$) and reweighting them by a factor of $\Pi^{-1}$ and $\Pi$, respectively. In fact, eq. 9 coincides with the Riccati equation associated with a slightly modified LQR problem:

$$ $\displaystyle\;\;\min_{u\in L^{2}[0,\infty)}\;\;\lim_{T\to\infty}\mathbb{E} \left[\frac{1}{T}\int_{0}^{T}x^{\top}\Pi^{-1}BB^{\top}\Pi^{-1}x+u^{\top}\Pi^{2 }u\;dt\right]$ (10) $$

This is the regulator problem for the adjoint state $x_{a}(t)=\Pi^{-1}x(t)$. Therefore, under the assumption that the observed dynamics can be approximated by a linear dynamical system,
we can measure LQG singular values associated with this modified LQR problem directly from measuring the causal minimum mean square error (MMSE) associated with prediction of $x(t)$ ($Q$), and the acausal MMSE associated with prediction of $x_{a}(t)$ ($\tilde{P}$).

To explicitly construct an estimator of the quantity $\text{Tr}(\tilde{P}Q)=\text{Tr}(Q\tilde{P})$, recall the matrix $Q$ is the error covariance of MMSE prediction of the system state $x(t)$ given past observations $y(t)$ over the interval $(t-T,t)$, whereas the matrix $\tilde{P}$ is the error covariance of MMSE prediction of the transformed system state $x_{a}(t)$ given future observations $y_{a}(t)$ over the interval $(t,t+T)$. The choice of $T$ is the only hyperparameter associated with FCCA. As discussed above, the Kalman Filter is used to efficiently calculate these MMSE estimates given an explicit state space model of the dynamics. In our case, to keep system dynamics implicit, we instead directly use the formulas for the MMSE error covariance in terms of cross correlations between $x(t),x_{a}(t)$ and $y(t),y_{a}(t)$. The standard formulas for the error covariance of MMSE prediction of a Gaussian distributed variable $z$ given $v$ read: $\Sigma_{z}-\Sigma_{zv}\Sigma_{v}^{-1}\Sigma_{vz}^{\top}$ where $\Sigma_{z}=\mathbb{E}[zz^{\top}],\Sigma_{v}=\mathbb{E}[vv^{\top}]$ and $\Sigma_{zv}=\mathbb{E}[zv^{\top}]$. The FCCA objective function is thus:

$$ $\displaystyle\text{FCCA}:\quad\text{argmin}_{C}\text{Tr}\left[\underbrace{\Big {(}\Pi-\Lambda_{1:T}(C)\Sigma^{-1}_{T}(C)\Lambda^{\top}_{1:T}(C)\Big{)}}_{ \text{causal MMSE covariance }(Q)}\underbrace{\Big{(}\Pi^{-1}-\tilde{\Lambda}_ {1:T}^{\top}(C)\Sigma^{-1}_{T}(C)\tilde{\Lambda}_{1:T}(C)\Big{)}}_{\text{ acausal MMSE covariance }(\tilde{P})}\right]$ (11) $$

where for discretization timescale $\tau$,

$$ $\displaystyle\Pi$ $\displaystyle=\underset{\text{(covariance of the neural data)}}{\mathbb{E}[x(t )x(t)^{\top}]},\;\Lambda_{k}=\underset{(\text{autocorrelation of the neural data})}{\mathbb{E}[x(t+k\tau)x(t)^{\top}]},\;\underset{(\text{autocorrelations of the adjoint state})}{\tilde{\Lambda}_{k}=\mathbb{E}[x_{a}(t+k\tau)x_{a}(t) ^{\top}]}$ $\displaystyle\Lambda_{1:T}(C)$ $\displaystyle=\{\Lambda_{1}C^{\top},\Lambda_{2}C^{\top},...,\Lambda_{T}C^{\top }\},\;\tilde{\Lambda}_{1:T}(C)=\{\tilde{\Lambda}_{1}\Pi C^{\top},\tilde{ \Lambda}_{2}\Pi C^{\top},...,\tilde{\Lambda}_{T}\Pi C^{\top}\}$ $$

and $\Sigma_{T}(C)$ is a block-Toeplitz space by time covariance matrix of $y(t)$ (i.e. the $ij^{\text{th}}$ block of $\Sigma_{T}(C)$ is given by $C^{\top}\Lambda_{|i-j|}C$. We optimize the FCCA objective function via L-BFGS.

### 2.4 Control-Theoretic Intuition for FCCA

We have shown how the sum of LQG singular values is an intrinsic measure of the cost to filter/regulate a linear dynamical system which is minimized at a fixed readout dimensionality by FCCA. We now provide further intuition for FCCA. In order to control the system state and carry out the computations necessary to perform state estimation and control signal synthesis, the controller itself must implement its own internal state dynamics. Thus, in addition to the complexity of the system itself, we may inquire about the complexity of the controller. One intuitive measure of this complexity is given by the controller’s state dimension (i.e., the McMillan degree), or the number of dynamical degrees of freedom it must implement to function. In the context of brain circuits, the degrees of freedom of the controller must ultimately be implemented via networks of neurons. We therefore hypothesize that biology favors performing task relevant computations via dynamics that require low dimensional controllers to regulate. As we argue below, minimizing the sum of LQG singular values over readout matrices ($C$) corresponds to a relaxation of the objective of searching for a subspace that enables control via a controller of low dimension. In other words, FCCA searches for dynamics that can be regulated with controllers of low complexity.

Figure: Figure 1: In principle, a controller of dimension as large as the neural state space may be required to effectively regulate dynamics within a FBC subspace ($H_{1}(s)$). However, subspaces optimized to minimize either the rank, or more practically, the trace of $PQ$ will require controllers of lower dimensionality to achieve near-optimal performance ($H_{2}(s)$).

Recall from Theorem 1 above that there exists a linear transformation that simultaneously diagonalizes both $P$ and $Q$. Let $(\tilde{A},\tilde{B},\tilde{C})$ be the corresponding
balanced realization. Order the LQG singular values in descending magnitude $\{\mu_{1},...,\mu_{N}\}$ and divide them into two sets $\{\mu_{1},...,\mu_{m}\}$ and $\{\mu_{m+1},...,\mu_{N}\}$. Assume the system input is of dimensionality $p$ and the output is of dimension $d$ (i.e., $\tilde{B}\in\mathbb{R}^{N\times p}$ and $\tilde{C}\in\mathbb{R}^{d\times N}$). Then, one can partition the state matrices $\{\tilde{A},\tilde{B},\tilde{C}\}$ accordingly:

$$ $\displaystyle\tilde{A}$ $\displaystyle=\begin{bmatrix}A_{11}&A_{12}\\ A_{21}&A_{22}\end{bmatrix}\quad\tilde{B}=\begin{bmatrix}B_{1}\\ B_{2}\end{bmatrix}\quad\tilde{C}=\begin{bmatrix}C_{1}&C_{2}\end{bmatrix}$ $$

Where $A_{11}\in\mathbb{R}^{m\times m},A_{22}\in\mathbb{R}^{N-m\times N-m}$, $B_{1}\in\mathbb{R}^{m\times p},B_{2}\in\mathbb{R}^{N-m\times p}$, $C_{1}\in\mathbb{R}^{d\times m},C_{2}\in\mathbb{R}^{d\times N-m}$. It can be shown that the optimal controller of dimension $m$ is obtained
from solving the Riccati equations corresponding to the truncated system $(A_{11},B_{1},C_{1})$. If the LQG singular values $\{\mu_{m+1},...,\mu_{N}\}$ are negligible, then the controller dimension can be reduced with essentially no loss in regulation performance. We illustrate this idea schematically in Figure 1, where the controller with transfer function $H_{1}(s)$ is approximated by a controller with lower state dimension $H_{2}(s)$. This suggests that to search for subspaces of neural dynamics that require low dimensional controllers to regulate, one should minimize the objective function $\text{argmin}_{C}\text{Rank}(\tilde{P}Q)$, where $\tilde{P}$ and $Q$ are the solutions to the Riccati equations 9 and 5, respectively. However, rank minimization is an NP-hard problem. A convex relaxation of the rank function is the nuclear norm (i.e. the sum of the singular values) [Fazel et al., ]. Given that $\tilde{P}Q$ is a positive semi-definite matrix, a tractable objective function that seeks subspaces of dynamics that require low complexity controllers is given by:

$$ $\displaystyle\text{argmin}_{C}\text{Tr}(\tilde{P}Q)$ $$

which is precisely what FCCA minimizes in a data-driven fashion (eq. 11).

## 3 PCA and FCCA subspaces diverge in non-normal dynamical systems

Figure: Figure 2: (Black) Average subspace angles between $d=2$ FCCA and PCA projections applied to Dale’s law constrained linear dynamical systems (LDS) as a function of non-normality. Spread indicates standard deviation over 20 random generations of A and 10 random initializations of FCCA. (Blue) Subspace angles between $d=2$ FCCA and PCA projections applied to firing rates derived from spiking activity driven by Dale’s Law constrained LDS. Spread around both curves indicates standard deviation taken over 20 random generations of $A$ matrices and 10 random initializations of FCCA.
Refer to caption: x2.png

Having derived data driven optimization problems to identify feedforward (PCA) and feedback (FCCA) controllable subspaces, we investigated under what conditions the solutions of PCA and FCCA will be distinct. We found that a key feature of the dynamical system of eq. 1 that determines the similarity of PCA and FCCA solutions is the non-normality of the underlying dynamics matrix, $A$. We first prove that when $A$ is normal (symmetric), and $B=I$, the critical points of PCA (eq. 11) and the FCCA objective function (eq. 7) coincide. (^1^11The set of real-valued, normal $A$ matrices can be divided into symmetric and orthogonal matrices. We restrict our treatment to stable dynamical systems. As orthogonal matrices give rise to systems that are only marginally stable, below we will use normal $A$ to refer interchangeably to symmetric $A$.)

###### Theorem 2 .

For $B=I_{N},A=A^{\top},A\in\mathbb{R}^{N\times N}$, with all eigenvalues of $A$ distinct and $\max\text{Re}(\lambda(A))<0$, the critical points of the feedforward controllability objective function eq. 4 and the feedback controllability objective function eq. 7 for projection dimension $d$ coincides with the eigenspace spanned by the $d$ eigenvalues with largest real value.

The proof of the theorem is provided in the Appendix. The restriction to $B=I$ is made within the proof, but does not apply to the general application of the method. Intuitively, in the case of symmetric, stable $A$, perturbations exponentially decay in all directions, and so the maximum response variance, and hence greatest feedforward controllability, is contained in the subspace with slowest decay, which corresponds to the eigenspace spanned by the $d$ eigenvalues with largest real value. The intuition for the slow eigenspace of $A$ serving as a (locally) optimal projection in the feedback controllability case is given by the fact that state reconstruction from past observations, the goal of the Kalman filter, will occur optimally using observations that have maximal autocorrelations with future state dynamics. Similarly, for the LQR, for a fixed rank input, the most variance will be suppressed by regulating within the subspace with slowest relaxation dynamics.

However, due to Dale’s Law, brain dynamics will be generated by non-normal dynamical systems. To demonstrate the effect of increasing the non-normality of $A$ on the solutions of PCA and FCCA, we turn to numerical simulations (the optimal feedback controllable projections are not analytically tractable). We generated 200-dimensional dynamics matrices constrained to follow Dale’s Law with an equal number of excitatory and inhibitory neurons. Neurons were connected randomly with a uniform connection probability of 0.25. To tune the non-normality of the system, we vary the strength of synaptic weights in the neuronal connectivity matrix. The strength of synaptic weights determines the spectral radius of the corresponding matrices [Rajan and Abbott, ]. Leaving the excitatory weights fixed, we then optimize the inhibitory weights as detailed in [Hennequin et al., ] to ensure system stability. The resulting matrices will have enhanced non-normality, with the degree of resulting non-normality having, empirically, a monotonic relationship with the starting spectral radius. We applied our methods both directly to the cross-covariance matrices of the resulting linear dynamical systems, as well as to spiking activity driven by simulated $x_{t}$. In the latter case, spiking activity was generated as a Poisson process with rate $\lambda_{t}=\exp(x_{t})$. Firing rates were obtained by binning spikes and applying a Gaussianizing boxcox transformation [Sakia, ]. These rates were then used to estimate the cross-covariance matrices. This procedure mirrors that which was applied to neural data in the subsequent section.

In Figure 2, we plot the average subspace angles between FCCA and PCA for $d=2$ projections (other choices of $d$ shown in Supplementary Figure 1 applied both directly to cross-covariance matrices of the linear dynamical systems (LDS, black) and cross-covariance matrices estimated from spiking activity (Count LDS, blue) as a function of the non-normality of the underlying $A$ matrix (measured using the Henrici metric, $||A^{\top}A-AA^{\top}||_{F}$). In both cases, we observe a nearly monotonic increase in the angles between FCCA and PCA subspaces as non-normality is increased. We note that as we constrain $A$ matrices to follow Dale’s Law, we cannot tune them to be completely normal, and hence the subspace angles between FCCA and PCA remain bounded away from zero even at the lower end of non-normality. In Supplementary Figure 2, we show that large subspace angles between FCCA and PCA also persist in the case when systems dynamics are made non-stationary by switching between a sequence of different linear systems. Thus, this new control-theoretic result suggests that PCA and FCCA subspaces should be geometrically distinct in neural population data, given the generality of non-normality of dynamics due to Dale’s Law.

###### Theorem 2 .

## 4 FCCA subspaces are better predictors of behavior than PCA subspaces

Figure: Figure 3: (a) Average subspace angles between FCCA and PCA at $d=2$ across recording sessions (median $\pm IQR$ indicated). (b) Five-fold cross-validated position prediction $r^{2}$ as a function of projection dimension between for FCCA (red) and PCA (black) and without dimensionality reduction (dashed blue). Mean $\pm$ standard error across folds and recording sessions indicated. (inset) Total area under the curve (AUC) of decoding performance averaged over folds for PCA and FCCA within each recording session (** : $p<10^{-2},n=8$, Wilcoxon signed rank test)

**Table 1: FCCA/PCA comparison across neural datasets**
| Dataset/Brain Region | $N_{r}$ | $\theta$ (deg.) | Peak Percent $\Delta$-$r^{2}$ | $\Delta$-$r^{2}$ AUC |
| --- | --- | --- | --- | --- |
| Hippocampus | 8 | $74.6\pm 1.7$ | $465\pm 144\%$ | $3.14\pm 0.30$ |
| M1 random | 35 | $58.0\pm 1.1$ | $229\pm 58\%$ | $2.75\pm 0.12$ |
| S1 random | 8 | $67.5\pm 3.8$ | $761\pm 189\%$ | $2.47\pm 0.36$ |
| M1 maze | 5 | $49.4\pm 4.3$ | $290\pm 72\%$ | $1.45\pm 0.23$ |

We first applied FCCA to tetrode recordings from the rat hippocampus made during a maze navigation
task. Further details on the dataset and preprocessing steps used are provided in the Appendix. In each recording session, we fit PCA and FCCA to neural activity across a range of projection dimensions. In line with the predictions of our theory and numerical simulations, we find that the subspace angle between PCA and FCCA was consistently large across recording sessions ($>3\pi/8$, Figure 3a, median and IQR indicated). We used $T=3$ (time bins) as the FCCA hyperparameter. Thus, we find that feedforward and feedback controllable subspaces are geometrically distinct in neural activity.

We next assessed the extent to which feedback controllable dynamics, as identified by FCCA, as opposed to feedforward controllable dynamics, as identified by PCA, were relevant for behavior. We trained linear decoders of the rat position from activity projected into FCCA and PCA subspaces. We used a window of 300 ms of neural activity centered around each time point to predict the corresponding binned position variable. We used linear decoders to emphasize the structure in the different subspaces available to a simple read-out. In Figure 3b, we report five-fold cross-validated prediction accuracy for PCA (black) and FCCA (red) over a range of projection dimensions (mean $\pm$ standard error across recording sessions and folds indicated). We found activity within FCCA subspaces to be more predictive of behavior than PCA subspaces across all dimensions, with a peak improvement of 112% at $d=13$. This superior decoding performance additionally held consistently across each recording session individually. In the inset of Fig. 3c, we plot the total area under prediction $r^{2}$ curves shown for each recording session (FCCA significantly higher than PCA, **: $p<10^{-2},n=8$, Wilcoxon signed rank test). FCCA is a nonconvex optimization problem. In practice, we optimize the FCCA objective functions over many randomly initialized orthogonal projection matrices and choose the final solution that yields the lowest value of the cost function 11. Feedback controllable subspaces therefore better capture behaviorally relevant dynamics than feedforward controllable subspaces.

To validate the robustness of these results, we repeated our analyses in two other datasets: recordings from macaque primary motor (M1 random) and primary somatosensory (S1 random) cortices during a self paced reaching task (O’Doherty et al. []), and recordings from macaque primary motor cortex during a delayed reaching task (M1 maze, Churchland et al. []). Further details on data preprocessing are provided in the Appendix. In Table 1, we report the number of recording sessions ($N_{r}$), mean: average subspace angle between FCCA and PCA subspaces at $d=2$ ($\theta$), peak percent $\Delta$-$r^{2}$ of behavioral prediction, and difference in the area under the behavioral prediction curves between PCA and FCCA. In all cases, standard errors are taken across the recording sessions, and analogously to Figure 3, behavioral decoding was performed from $d=1$ to $d=30$. Importantly, in all datasets, FCCA performed better behavioral prediction, and the subspace angles between FCCA and PCA were substantially different from zero. In Supplementary Figure 3, we confirm that the substantial subspace angles between FCCA and PCA are largely insensitive to the choice of T, the choice of projection dimensionality, and robust across initializations of FCCA. Further, in Supplementary Figure 4, we verify that the superior decoding performance of FCCA subspaces hold consistently across each individual initialization.

## 5 Discussion

We presented FCCA, a novel linear dimensionality reduction method that identifies feedback controllable subspaces of neural population dynamics. The correspondence between PCA and feedforward controllability, long known in the control theory community [Moore, ], but unrecognized in the neuroscience community, adds additional interpretative value to these subspaces. We demonstrated that feedforward and feedback controllable subspaces are geometrically distinct in non-normal dynamical systems, a fact of fundamental importance to the analysis of neural dynamics from cortex, where Dale’s Law necessitates non-normality. Correspondingly, in electrophysiology recordings, we found large subspace angles between FCCA and PCA subspaces. Furthermore, we found that FCCA subspaces were better predictors of behavior than PCA subspaces. This suggests that targeting feedback controllable subspaces in the design of brain machine interfaces may be advantageous in terms of accuracy of behavioral prediction, the number of samples needed to calibrate predictions to a desired level of accuracy, and the efficacy of closed loop perturbations.

Several methodological extensions to FCCA are possible. In FCCA, we rely on estimation of the regulator cost through acausal filtering (eq. 9 and estimate the filtering error through the Gaussian MMSE formula (eq. 11) to keep the model of the data implicit. These correspondences only hold for linear systems under a particular choice of the LQR cost function (eq. 10). While this makes the method computationally efficient, it restricts the form of weight matrices in the LQR objective functions that can be considered. The objective function in eq. 7 could alternatively be applied to post-hoc analysis of linear state space models fit to neural recordings [Gao et al., ], as these models explicitly yield the system matrices required to solve the Riccati equations 5 and 6. This analysis could be combined with techniques from inverse linear optimal control [Priess et al., ] to provide a more refined picture of the controllability of population dynamics. Generalization to nonlinear measures of controllability Scherpen [], Bouvrie and Hamzi [] is another direction of future work.

## Appendix A Appendix

### A.1 Details of neural datasets

Data from the hippocampus contained recordings from a single rodent. There were a total of 8 recording sessions lasting approximately 20 minutes each with between 98-120 identified single units within each recording session. We performed our analyses on neural activity while the rat was in motion (velocity > 4 cm/s).

The M1/S1 random dataset contained a total of 35 recording sessions from 2 monkeys (28 within monkey 1, 7 within monkey 2) spanning 17309 total reaches (13149 from monkey 1, 4160 from monkey 2). Of the 35 recording sessions, 8 included activity from S1. The number of single units in each recording session varied between 96-200 units in M1, and 86-187 in S1. The maze dataset contained 5 recording sessions recorded from 2 different monkeys comprising 10829 total reaches (8682 in monkey 3, 2147 in monkey 4). Each recording session contained 96 single units. Both datasets mapped the monkey hand location to a cursor location on the 2D task plane. For the M1/S1 random dataset, we decoded cursor velocity, whereas for the maze dataset, we decoded cursor position.

We binned spikes within the hippocampal data at 25 ms, and the M1/S1 random and M1 maze datasets at 50 ms. We then applied a boxcox transformation to binned firing rates to Gaussianize the data. A single fit of FCCA on the activity from a single recording session in the datasets considered using a desktop computer equipped with an 8 core CPU and 64 GB of memory requires < 5 seconds.

### A.2 Proof of Theorem 2

In this section, we prove the equivalence of the solutions of the FFC (eq. 4) and FBC objective functions (eq. 7) when system dynamics are stable and symmetric. We focus on symmetric matrices as the requirement that dynamics be stable (i.e., all eigenvalues of the dynamics $A$ must have negative real part) essentially reduces the space of normal matrices to that of symmetric matrices. We reproduce these objective functions for convenience:

$$ $\displaystyle C_{\text{FFC}}:\;\;\text{argmax}_{C}\log\det C\Pi C^{\top}$ $\displaystyle C_{\text{FBC}}:\;\;\text{argmin}_{C}\text{Tr}(PQ)$ $$

We prove this theorem when the matrix $P$ in the FBC objective function arises from the canonical LQR loss function:

$$ $\min_{u}\left\{\lim_{T\to\infty}\mathbb{E}\left[\frac{1}{T}\int_{0}^{T}x^{\top }x+u^{\top}u\;dt\right],\;\;x(0)=x_{0},u\in L^{2}[0,\infty)\right\}$ $$

and not the variant given in eq. 10. When calculating FBC from data within FCCA, we must use the latter LQR loss function as it maps onto acausal filtering, and therefore may be estimated from data. Recall from the discussion below eq. 8 that within the FFC objective function, we assess controllability when the output/observation matrix $C$ is used as the input matrix for the regulator signal (i.e., we make the relabeling $B^{\top}\to C$. We further work under the assumption that the input matrix $B$ to the open loop system is equal to the identity. The open loop dynamics of $x(t)$ are then given by:

$$ $\displaystyle\dot{x}=Ax(t)+u(t)$ (12) $$

where $u(t)$ has the same dimensionality as $x(t)$, and is uncorrelated with the past of $x(t)$ (i.e. $u(t)\perp x(\tau),\tau<t$). Formally, $u(t)$ represents the innovations process of $x(t)$. The equations for $Q$ (corresponding to the Kalman Filter, eq. 5) and the equation for $P$ (corresponding to the LQR, eq. 6) reduce to the following:

$$ $\displaystyle AQ+QA+I_{N}-QC^{\top}CQ$ $\displaystyle=0$ (13) $\displaystyle AP+PA+I_{N}-PCC^{\top}P$ $\displaystyle=0$ (14) $$

where $I_{N}$ denotes the $N\times N$ identity matrix.

We observe that under the stated assumptions, the Riccati equations for $Q$ and $P$ actually coincide, and thus the FBC objective function reads $Tr(Q^{2})$. We will show that both FFC and FBC objective functions achieve local optima for some fixed projection dimension $d$ when the projection matrix $C$ coincides with a projection onto the eigenspace spanned by the $d$ eigenvalues of A with largest real part, which we denote as $V_{d}$. In fact, in the case of the FFC objective function, the eigenspace corresponds to a global optimum. For the FBC objective function, we are able to establish global optimality rigorously only for the $2D\to 1D$ dimension reduction.

We briefly outline the proof strategy. First, we will prove the optimality of $V_{d}$ for the FFC objective function in section S1.9.1 by showing that (i) $V_{d}$ is an eigenvector of $\Pi$ in the case when $A$ is symmetric and (ii) relying on the Poincare Separation Theorem. Then, in section S1.9.2, we will prove that $V_{d}$ is a critical point of the FBC objective function. The proof relies on an iterative technique to solve the Riccati equation. These iterates form a recursively defined sequence that provide increasingly more accurate approximations to the FBC objective function that converge in the limit. Treating these iterative approximations of the FBC objective function as a function of $C$, we show that $V_{d}$ is a critical point of all iterates, and thus in the limit, $V_{d}$ is a critical point of the FBC objective function.

#### FFC Objective Function

###### Lemma 1 .

For $B=I_{N},A=A^{\top},A\in\mathbb{R}^{N\times N}$, with all eigenvalues of $A$ distinct and $\max\text{Re}(\lambda(A))<0$, the optimal solution for the feedforward controllability objective function for projection dimension $d$ coincides with $V_{d}$, the matrix whose rows are formed by the eigenvectors corresponding to the $d$ eigenvalues of $A$ with largest real value.

*Proof*

The FFC objective function reads:

$$ $\displaystyle\text{argmax}_{C}\log\det C\Pi C^{\top}\;\;|\;\;C\in\mathbb{R}^{d \times N},CC^{\top}=I_{d}$ (15) $$

We first re-write $\Pi$:

$$ $\displaystyle\Pi=\int_{0}^{\infty}dt\;e^{At}BB^{\top}e^{A^{\top}t}=\int_{0}^{ \infty}dt\;e^{2At}$ $$

Let $A=U\Lambda U^{\top}$ denote the eigenvalue decomposition of $A$. Recall that since $A=A^{\top}$, $U$ is orthogonal. Then we can write:

$$ $\displaystyle\Pi$ $\displaystyle=U\int_{0}^{\infty}dte^{2\Lambda t}U^{\top}$ $\displaystyle=\frac{1}{2}UDU^{\top}$ $$

where $D$ is a diagonal matrix with diagonal entries $\{\frac{1}{-\lambda_{1}},\frac{1}{-\lambda_{2}},...,\frac{1}{-\lambda_{N}}\}$ being the eigenvalues of $\Pi$. We conclude that the matrix $\Pi$ has the same eigenbasis as $A$. Also, since all $\lambda_{j}$ are real and negative, the ordering of the eigenvalues is preserved ($\lambda_{i}>\lambda_{j}$ implies $-\frac{1}{\lambda_{i}}>-\frac{1}{\lambda_{j}}$).

That $V_{d}$ solves 15 follows from the Poincare separation theorem, which we restate for convenience:

###### Proposition 1 .

Poincare Separation Theorem (Magnus and Neudecker [], 11.10)

Let $M$ be any square, symmetric matrix, and let $\mu_{1}\geq\mu_{2}\geq\ldots\geq\mu_{N}$ be its eigenvalues. Let $C\in\mathbb{R}^{d\times N}$ be a semi-orthogonal matrix (i.e., $CC^{\top}=I_{d}$). Then, the eigenvalues $\eta_{1}\geq\eta_{2}\geq\ldots\geq\eta_{d}$ of $CMC^{\top}$ satisfy:

$$ $\mu_{i}\geq\eta_{i}\geq\mu_{N-d+i}$ $$

In particular, Proposition 1 implies that $\det CMC^{\top}=\prod_{i=1}^{d}\eta_{i}\leq\prod_{i=1}^{d}\mu_{i}$, and hence $\log\det CMC^{\top}\leq\sum_{i=1}^{d}\log\mu_{i}$ We now show that this inequality is satisfied with equality when $C=V_{d}$. Consider the optimization problem

$$ $\displaystyle\text{argmax}_{C}\log\det CMC^{\top}\;\;|\;\;C\in\mathbb{R}^{d \times N},CC^{\top}=I_{d}$ (16) $$

Let $M=U\Gamma U^{\top}$ be the eigendecomposition of $M$. We can equivalently parameterize the optimization problem as:

$$ $\displaystyle\text{argmax}_{\tilde{C}}\log\det\tilde{C}\Gamma\tilde{C}^{\top} \;\;|\;\;\tilde{C}\in\mathbb{R}^{d\times N},\tilde{C}\tilde{C}^{\top}=I_{d}$ (17) $$

The solution to the original problem, eq. 16, can be recovered from setting $C=\tilde{C}U^{\top}$. Now, assume (without loss of generality) that we have arranged the values of $\Gamma$ so that the largest d eigenvalues, $\mu_{1},...,\mu_{d}$, occur first. We observe that the choice of $\tilde{C}=[I_{d};\mathbf{0}_{N-d,N-d}]\equiv\tilde{C}_{*}$, which picks out these first d elements of the diagonal of $\Gamma$, yields $\log\det\tilde{C}^{\top}_{*}\Gamma\tilde{C}_{*}=\sum_{i=1}^{d}\log\mu_{i}$, and hence solves the desired optimization problem. It follows that $C_{*}=\tilde{C}_{*}U^{\top}=V_{d}$

To complete the proof of Lemma 1, we substitute $M$ with $\Pi$, and the eigenvalues $\mu_{i}$ with $-1/\lambda_{i}$ (the eigenvalues of $\Pi$, expressed in terms of the eigenvalues of $A$).

$\;\square$
FBC Objective Function
For the case of the FBC objective function, we show that projection matrices of rank $d$ that align with the $d$ slowest eigenmodes of $A$ constitute local minima of the objective function. We rely on two simplifying features of the problem. First, the FBC objective function is invariant to the choice of basis in the state space. We therefore work within the eigenbasis of $A$, as within this basis, the system defined by eq. 12 decouples into $n$ non-interacting scalar dynamical systems. Additionally, we rely on the fact that the FBC objective function is also invariant to coordinate transformations within the projected space. In other words, the choice of coordinates in which we express $y$ also makes no difference. Without loss of generality then, we may treat the problem in a basis where $A$ is diagonal with entries given by its eigenvalues and $C$ is an orthonormal projection matrix (i.e. $CC^{\top}=I_{d}$). A restatement of the latter condition is that $C$ belongs to the Steifel manifold of $N\times d$ matrices: $\Omega\equiv\{C\in\mathbb{R}^{N\times d}|CC^{\top}=I_{d}\}$.

###### Lemma 2 .

For $B=I_{N},A=A^{\top},A^{N\times N}$, with all eigenvalues of $A$ distinct and $\max\text{Re}(\lambda(A))<0$, the projection matrix onto the eigenspace spanned by the $d$ eigenvalues of $A$ with largest real value constitutes a critical point of the LQG trace objective function on $\Omega$

*Proof*
Explicitly calculating the gradient of the solution of the Riccati equation is analytically intractable for $n>1$, and so we we will rely on the analysis of an iterative procedure to solve the Riccati equation via Newton’s method, known as the Newton-Kleinmann (NK) iterations [Kleinman, ]. These iterations are described in the following proposition:

###### Proposition 2 .

Consider the Riccati equation $0=AQ+QA^{\top}+BB^{\top}-QC^{\top}CQ$. Let $Q_{m},m=1,2,...$ be the unique positive definite solution of the Lyapunov equation:

$$ $\displaystyle 0=A_{k}Q_{m}+Q_{m}A_{k}^{\top}+BB^{\top}+Q_{m-1}C^{\top}CQ_{m-1}$ (18) $$

where $A_{k}=A-C^{\top}CQ_{m-1}$, and where $Q_{0}$ is chosen such that $A_{1}$ is a stable matrix (i.e. all real parts of its eigenvalues are $<0$). For two positive semidefinite matrices $M,N$, we denote $M\geq N$ if the difference $M-N$ remains positie semidefinite. Then:

- 1.
$Q\leq Q_{m+1}\leq Q_{m}\leq...,k=0,1$
- 2.
$\lim_{k\to\infty}Q_{m}=Q$

Thus the $Q_{m}$ iteratively approach the solution of the Riccati equation from above. Since in our case, the Riccati equations for $P$ and $Q$ coincide, an identical sequence $P_{k}$ can be constructed using analogous NK iterations that approaches $P$ from above. From this, it follows that $\lim_{k\to\infty}Tr(Q_{m}P_{k})=\lim_{k\to\infty}Tr(Q_{m}^{2})=Tr(Q^{2})$. We then use the fact that in addition to the $Q_{m}$ converging to $Q$, the sequence $\nabla_{C}\text{Tr}\left(Q_{m}^{2}\right)$ converges to $\nabla_{C}\text{Tr}(Q^{2})$ as $k\to\infty$, where $\nabla_{C}$ denotes the gradient with respect to $C$. This is rigorously established in the following lemma, which is the multivariate generalization of Theorem 7.17 from [Rudin and others, ]:

###### Lemma 3 .

Suppose $\{f_{m}\}$ is a sequence of functions differentiable on an interval $h\subset\mathcal{H}$, where $\mathcal{H}$ is some finite-dimensional vector space, such that $\{f_{m}(x_{0})\}$ converges for some point $x_{0}\in h$. If $\{\nabla f_{m}(x_{0})\}$ converges uniformly in $h$, then $\{f_{m}\}$ converges uniformly on $I$, to a function $f$, and

$$ $\nabla f(x)=\lim_{m\to\infty}\nabla f_{m}(x)\quad x\in h$ $$

Here, the $\{f_{m}\}$ are the Newton-Kleinmann iterates $Q_{m}$, and $x_{0}$ corresponds to the $C$ matrix that projects onto the slow eigenspace of $A$. The NK iterates are known to converge uniformly over an interval of possible $C$ matrices (in fact any such $C$ matrix for which there exists a $K$ such that $A-C^{T}CK$ is a stable matrix) [Kleinman, ].

We will calculate the gradient $\nabla_{C}Q_{m}$ on $\Omega$ by explicitly calculating the directional derivatives of $Q_{m}$ over a basis of the tangent space of $\Omega$ at $C_{\text{slow}}$. Any element $\Psi$ belonging to the tangent space at $C\in\Omega$ can be parameterized by the following [Edelman et al., ]:

$$ $\Psi=CM+(I_{N}-CC^{\top})T$ $$

where $M$ is skew symmetric and $t$ is arbitrary. Let $C_{\text{slow}}$ be the projection matrix onto the slow eigenspace of $A$ of dimension d. Since we work in the eigenbasis of $A$, $C_{\text{slow}}=\begin{bmatrix}I_{d}&0\end{bmatrix}$. At this point, elements of the tangent space take on the particularly simple form

$$ $\Psi=\begin{bmatrix}M&T\end{bmatrix}$ $$

where now $M$ is a $d\times d$ skew symmetric matrix and $T\in\mathbb{R}^{d\times(N-d)}$ is arbitrary. A basis for the tangent space is provided by the set of matrices $\{M_{ij},T_{kl},i=2,...d,j=1,...,i-1,k=1,...,d,l=1,...,N-d\}$ where $M_{ij}$ is a matrix with entry $1$ at index $(i,j)$ and $-1$ at index $(j,i)$ and zero otherwise, and $T_{kl}$ is the matrix with entry $1$ at index $(k,l)$ and zero otherwise. Denote by $D_{\Psi}Q_{m}$ the directional derivative of $Q_{m}$ along the direction of $\Psi$, viewing $Q_{m}$ as a function of C (denoted $Q_{m}[C]$):

$$ $\displaystyle D_{\Psi}Q_{m}=\lim_{\alpha\to 0}\frac{Q_{m}[C_{\text{slow}}+ \alpha\Psi]-Q_{m}[C_{\text{slow}}]}{\alpha}$ (19) $$

Let $\Psi_{ij,kl}$ denote the tangent matrix $\begin{bmatrix}M_{ij}&T_{kl}\end{bmatrix}$. Before calculating $Q_{m}(C_{\text{slow}}+\alpha\Psi_{ij,kl})$ explicitly, we first observe that as long as the NK iterations are initialized with a diagonal $Q_{0}$, then the diagonal nature of $C_{\text{slow}}^{\top}C_{\text{slow}}$ ensures that all $Q_{m}$ will subsequently remain diagonal matrices. In fact, it can be shown that $\lim_{k\to\infty}Q_{m}=Q$ will also be diagonal, in this case. We write $A$ in block form as $\begin{bmatrix}\Lambda_{||}&0\\
0&\Lambda_{\perp}\end{bmatrix}$, and similarly $Q_{m-1}=\begin{bmatrix}\mathcal{Q}_{||}&0\\
0&\mathcal{Q}_{\perp}\end{bmatrix}$, where $\Lambda_{||},\mathcal{Q}_{||}$ are $d\times d$ diagonal matrices defined on the image of $C_{\text{slow}}$ and $\Lambda_{\perp},\mathcal{Q}_{\perp}$ are diagonal matrices defined on the kernel of $C_{\text{slow}}$. We denote the individual diagonal elements of $\Lambda_{||},\mathcal{Q}_{||}$ as $\lambda_{i},\mathcal{Q}_{i},i=1,...,d$ and of $\Lambda_{\perp},\mathcal{Q}_{\perp}$ as $\lambda_{i},\mathcal{Q}_{i},i=d,...,N-d$. Then, equation 18 becomes:

$$ $\displaystyle\left(\begin{bmatrix}\Lambda_{||}&0\\ 0&\Lambda_{\perp}\end{bmatrix}-\begin{bmatrix}(I_{d}-\alpha^{2}M_{ij}^{2}) \mathcal{Q}_{||}&(\alpha T_{kl}+\alpha^{2}M_{ij}^{\top}T_{kl})\mathcal{Q}_{ \perp}\\ (\alpha T_{kl}^{\top}+\alpha^{2}T_{kl}^{\top}M_{ij})\mathcal{Q}_{||}&\alpha^{2 }T_{kl}^{\top}T_{kl}\mathcal{Q}_{\perp}\end{bmatrix}\right)Q_{m}[C_{\text{slow }}+\Psi_{ij,kl}]$ $\displaystyle+Q_{m}[C_{\text{slow}}+\Psi_{ij,kl}]\left(\begin{bmatrix}\Lambda_ {||}&0\\ 0&\Lambda_{\perp}\end{bmatrix}-\begin{bmatrix}\mathcal{Q}_{||}(I_{d}-\alpha^{2 }M_{ij}^{2})&\mathcal{Q}_{||}(\alpha T_{kl}+\alpha^{2}M_{ij}^{\top}T_{kl})\\ \mathcal{Q}_{\perp}(\alpha T_{kl}^{\top}+\alpha^{2}T_{kl}^{\top}M_{ij})& \mathcal{Q}_{\perp}\alpha^{2}T_{kl}^{\top}T_{kl}\end{bmatrix}\right)$ $\displaystyle+I_{N}+\begin{bmatrix}\mathcal{Q}_{||}(I_{d}-\alpha^{2}M_{ij}^{2} )\mathcal{Q}_{||}&\mathcal{Q}_{||}(\alpha T_{kl}+\alpha^{2}M_{ij}^{\top}T_{kl} )\mathcal{Q}_{\perp}\\ \mathcal{Q}_{\perp}(\alpha T_{kl}^{\top}+\alpha^{2}T_{kl}^{\top}M_{ij}) \mathcal{Q}_{||}&\alpha^{2}\mathcal{Q}_{\perp}T_{kl}^{\top}T_{kl}\mathcal{Q}_{ \perp}\end{bmatrix}=0$ (20) $$

where we have used $M^{\top}=-M$. The equivalent equation for $Q_{m}(C_{\text{slow}})$ reads:

$$ $\displaystyle\left(\begin{bmatrix}\Lambda_{||}&0\\ 0&\Lambda_{\perp}\end{bmatrix}-\begin{bmatrix}\mathcal{Q}_{||}&0\\ 0&0\end{bmatrix}\right)Q_{m}[C_{\text{slow}}]+Q_{m}[C_{\text{slow}}]\left( \begin{bmatrix}\Lambda_{||}&0\\ 0&\Lambda_{\perp}\end{bmatrix}-\begin{bmatrix}\mathcal{Q}_{||}&0\\ 0&0\end{bmatrix}\right)+I_{N}+$ (21) $\displaystyle\begin{bmatrix}\mathcal{Q}_{||}^{2}&0\\ 0&0\end{bmatrix}=0$ (22) $$

This latter equation is easily solved to yield:

$$ $Q_{m}[C_{\text{slow}}]=\begin{bmatrix}\frac{1}{2}\left(I_{d}+\mathcal{Q}_{||}^ {2}\right)\left(\mathcal{Q}_{||}-\Lambda_{||}\right)^{-1}&0\\ 0&-\frac{1}{2}\Lambda_{\perp}^{-1}\end{bmatrix}$ $$

To explicitly solve the former equation, we recall that the matrices $M_{ij}$ and $T_{kl}$ have only two and one nonzero terms, respectively. $M_{ij}^{2}$ contains two nonzero terms at index $(i,i)$ and $(j,j)$. $T_{kl}^{\top}T_{kl}$ contains one non-zero term at index $(l,l)$. $M_{ij}^{\top}T_{kl}$ contains a single nonzero term at $(i,l)$ or $(j,l)$ only if $k=i$ or $k=j$, respectively. Accordingly, we distinguish between where $k=i$ or $k=j$ (without loss of generality we may assume that $k=j$), and where $k\neq i$ and $k\neq j$.

In what follows, we will denote the $(i,j)$ entry of $Q_{m}[C_{\text{slow}}+\alpha\Psi_{ij,kl}]$ as $q_{ij}$.

- 1.
*Case 1: $k=j$*
In this case, careful inspection of eq. 20 reveals that it differs from eq. 22 only within a $3\times 3$ subsystem:
$\displaystyle\begin{bmatrix}\mathcal{S}_{11}&\mathcal{S}_{12}&\mathcal{S}_{13}
\\
\mathcal{S}_{21}&\mathcal{S}_{22}&\mathcal{S}_{23}\\
\mathcal{S}_{31}&\mathcal{S}_{32}&\mathcal{S}_{33}\end{bmatrix}=0$
Note that this matrix is symmetric, yielding 6 equations for 6 unknowns:
$\displaystyle\mathcal{S}_{11}$
$\displaystyle=\alpha^{2}\mathcal{Q}_{i}^{2}+2\alpha^{2}\mathcal{Q}_{d+l}q_{i,d
+l}+\mathcal{Q}_{i}^{2}+2q_{ii}\left(-\alpha^{2}\mathcal{Q}_{i}+\lambda_{i}-
\mathcal{Q}_{i}\right)+1$
$\displaystyle\mathcal{S}_{12}$
$\displaystyle=\alpha^{2}\mathcal{Q}_{d+l}q_{j,d+l}-\alpha\mathcal{Q}_{d+l}q_{i
,d+l}+q_{ij}\left(-\alpha^{2}\mathcal{Q}_{i}+\lambda_{i}-\mathcal{Q}_{i}\right
)+q_{ij}\left(-\alpha^{2}\mathcal{Q}_{j}+\lambda_{j}-\mathcal{Q}_{j}\right)$
$\displaystyle\mathcal{S}_{13}$
$\displaystyle=-\alpha^{2}\mathcal{Q}_{i}\mathcal{Q}_{d+l}+\alpha^{2}\mathcal{Q
}_{i}q_{ii}+\alpha^{2}\mathcal{Q}_{d+l}q_{d+l}-\alpha\mathcal{Q}_{j}q_{ij}+$
$\displaystyle q_{i,d+l}\left(-\alpha^{2}\mathcal{Q}_{d+l}+\lambda_{d+l}\right)
+q_{i,d+l}\left(-\alpha^{2}\mathcal{Q}_{i}+\lambda_{i}-\mathcal{Q}_{i}\right)$
$\displaystyle\mathcal{S}_{22}$
$\displaystyle=\alpha^{2}\mathcal{Q}_{j}^{2}-2\alpha\mathcal{Q}_{d+l}q_{j,d+l}+
\mathcal{Q}_{j}^{2}+2q_{jj}\left(-\alpha^{2}\mathcal{Q}_{j}+\lambda_{j}-
\mathcal{Q}_{j}\right)+1$
$\displaystyle\mathcal{S}_{23}$
$\displaystyle=\alpha^{2}\mathcal{Q}_{i}q_{ij}+\alpha\mathcal{Q}_{j}\mathcal{Q}
_{d+l}-\alpha\mathcal{Q}_{j}q_{jj}-\alpha\mathcal{Q}_{d+l}q_{d+l,d+l}+q_{j,d+l
}\left(-\alpha^{2}\mathcal{Q}_{d+l}+\lambda_{d+l}\right)+$
$\displaystyle q_{j,d+l}\left(-\alpha^{2}\mathcal{Q}_{j}+\lambda_{j}-\mathcal{Q
}_{j}\right)$
$\displaystyle\mathcal{S}_{33}$
$\displaystyle=2\alpha^{2}\mathcal{Q}_{i}q_{i,d+l}+\alpha^{2}\mathcal{Q}_{d+l}^
{2}-2\alpha\mathcal{Q}_{j}q_{j,d+l}+2q_{d+l,d+l}\left(-\alpha^{2}\mathcal{Q}_{
d+l}+\lambda_{d+l}\right)+1$
Direct solution is still infeasible, but noting our interest is in the behavior of solutions as $\alpha\to 0$, and only terms of $O(\alpha)$ will survive in the limit in eq. 19, we consider solving these equations perturbatively. That is, we express each $q_{ij}$ in a power series in $\alpha$: $q_{ij}=q_{ij}^{(0)}+q_{ij}^{(1)}\alpha+O(\alpha^{2})$. One obtains each coefficient in the expansion by plugging this form into the above matrix and setting all terms of the corresponding order in $\alpha$ to 0. The lowest order term, $q_{ij}^{(0)}$, coincides with the solution of the unperturbed system, eq. 22. Plugging in the expansion into the $3\times 3$ subsystem above, as well as the solution of the unperturbed system, and collecting all coefficients proportional to $\alpha$ yields the following system of equations:
$\displaystyle\begin{bmatrix}\mathcal{S}^{(1)}_{11}&\mathcal{S}^{(1)}_{12}&
\mathcal{S}^{(1)}_{13}\\
\mathcal{S}^{(1)}_{21}&\mathcal{S}^{(1)}_{22}&\mathcal{S}^{(1)}_{23}\\
\mathcal{S}^{(1)}_{31}&\mathcal{S}^{(1)}_{32}&\mathcal{S}^{(1)}_{33}\end{
bmatrix}=0$
$\displaystyle\mathcal{S}^{(1)}_{11}$
$\displaystyle=2\lambda_{i}q_{ii}^{(1)}-2\mathcal{Q}_{i}q_{ii}^{(1)}$
$\displaystyle\mathcal{S}^{(1)}_{12}$
$\displaystyle=\lambda_{i}q_{ij}^{(1)}+\lambda_{j}q_{ij}^{(1)}-\mathcal{Q}_{i}q
_{ij}^{(1)}-\mathcal{Q}_{j}q_{ij}^{(1)}$
$\displaystyle\mathcal{S}^{(1)}_{13}$
$\displaystyle=\lambda_{i}q_{i,d+l}^{(1)}+\lambda_{d+l}q_{i,d+l}^{(1)}-\mathcal
{Q}_{i}q_{i,d+l}^{(1)}$
$\displaystyle\mathcal{S}^{(1)}_{22}$
$\displaystyle=2\lambda_{j}q_{jj}^{(1)}-2\mathcal{Q}_{j}q_{jj}^{(1)}$
$\displaystyle\mathcal{S}^{(1)}_{23}$
$\displaystyle=\lambda_{j}q_{j,d+l}^{(1)}+\lambda_{d+l}q_{j,d+l}^{(1)}+\mathcal
{Q}_{j}\mathcal{Q}_{d+l}-\mathcal{Q}_{j}q_{j,d+l}^{(1)}-\frac{\mathcal{Q}_{j}
\left(\mathcal{Q}_{j}^{2}+1\right)}{-2\lambda_{j}+2\mathcal{Q}_{j}}+\frac{
\mathcal{Q}_{d+l}}{2\lambda_{d+l}}$
$\displaystyle\mathcal{S}^{(1)}_{33}$
$\displaystyle=2\lambda_{d+l}q_{d+l}^{(1)}$
Solving this system yields the following solutions for the $q_{ij}^{(1)}$:
$\displaystyle q_{ii}^{(1)}$
$\displaystyle=0$
$\displaystyle q_{jj}^{(1)}$
$\displaystyle=0$
$\displaystyle q_{d+l,d+l}^{(1)}$
$\displaystyle=0$
$\displaystyle q_{ij}^{(1)}$
$\displaystyle=0$
$\displaystyle q_{i,d+l}^{(1)}$
$\displaystyle=0$
$\displaystyle q_{j,d+l}^{(1)}$
$\displaystyle=\frac{-2\lambda_{j}\lambda_{d+l}\mathcal{Q}_{j}\mathcal{Q}_{d+l}
-\lambda_{j}\mathcal{Q}_{d+l}-\lambda_{d+l}\mathcal{Q}_{j}^{3}+2\lambda_{d+l}
\mathcal{Q}_{j}^{2}\mathcal{Q}_{d+l}-\lambda_{d+l}\mathcal{Q}_{j}+\mathcal{Q}_
{j}\mathcal{Q}_{d+l}}{2\lambda_{j}^{2}\lambda_{d+l}+2\lambda_{j}\lambda_{d+l}^
{2}-4\lambda_{j}\lambda_{d+l}\mathcal{Q}_{j}-2\lambda_{d+l}^{2}\mathcal{Q}_{j}
+2\lambda_{d+l}\mathcal{Q}_{j}^{2}}$
- 2.
*Case 2*: $k\neq i,k\neq j$. In this case, we must again consider the $3\times 3$ subsystem indexed by $i,j,d+l$, but since $M_{ij}T_{kl}$ is a matrix of all zeros, the expression simplifies considerably:
$\displaystyle\begin{bmatrix}\mathcal{S}_{11}&\mathcal{S}_{12}&\mathcal{S}_{13}
\\
\mathcal{S}_{21}&\mathcal{S}_{22}&\mathcal{S}_{23}\\
\mathcal{S}_{31}&\mathcal{S}_{32}&\mathcal{S}_{33}\end{bmatrix}=0$
$\displaystyle\mathcal{S}_{11}$
$\displaystyle=\alpha^{2}\mathcal{Q}_{i}^{2}+\mathcal{Q}_{i}^{2}+2q_{i}\left(-
\alpha^{2}\mathcal{Q}_{i}+\lambda_{i}-\mathcal{Q}_{i}\right)+1$
$\displaystyle\mathcal{S}_{12}$
$\displaystyle=q_{ij}\left(-\alpha^{2}\mathcal{Q}_{i}+\lambda_{i}-\mathcal{Q}_{
i}\right)+q_{ij}\left(-\alpha^{2}\mathcal{Q}_{j}+\lambda_{j}-\mathcal{Q}_{j}\right)$
$\displaystyle\mathcal{S}_{13}$
$\displaystyle=\lambda_{d+l}q_{i,d+l}+q_{i,d+l}\left(-\alpha^{2}\mathcal{Q}_{i}
+\lambda_{i}-\mathcal{Q}_{i}\right)$
$\displaystyle\mathcal{S}_{22}$
$\displaystyle\alpha^{2}\mathcal{Q}_{j}^{2}+\mathcal{Q}_{j}^{2}+2q_{j}\left(-
\alpha^{2}\mathcal{Q}_{j}+\lambda_{j}-\mathcal{Q}_{j}\right)+1$
$\displaystyle\mathcal{S}_{23}$
$\displaystyle=\lambda_{d+l}q_{j,d+l}+q_{j,d+l}\left(-\alpha^{2}\mathcal{Q}_{j}
+\lambda_{j}-\mathcal{Q}_{j}\right)$
$\displaystyle\mathcal{S}_{33}$
$\displaystyle=2\lambda_{d+l}q_{d+l}+1$
Plugging in the power series expansion $q_{ij}=q_{ij}^{(0)}+q_{ij}^{(1)}\alpha+O(\alpha^{2})$, one finds the lowest order terms in $\alpha$ within this system of equations occurs at $O(\alpha^{2})$, and thus to $O(\alpha)$, the solution of $Q_{m}[C_{\text{slow}}+\alpha\Psi_{ij,kl}]$ coincides with $Q_{m}[C_{\text{slow}}]$.

To complete the proof of Theorem 3, we must calculate the following quantity:

$$ $\displaystyle D_{\Psi_{ij,kl}}\text{Tr}\left(Q_{m}^{2}\right)=\lim_{\alpha\to 0 }\frac{\text{Tr}(Q_{m}[C_{\text{slow}+\alpha\Psi_{ij,kl}}]^{2})-\text{Tr}(Q_{m }[C_{\text{slow}}]^{2})}{\alpha}$ $$

From the case-wise analysis above, we see that the only matrix element of $Q_{m}$ that differs between $Q_{m}[C_{\text{slow}+\alpha\Psi_{ij,kl}}]$ and $Q_{m}[C_{\text{slow}}]$ to $O(\alpha)$ is an off-diagonal term ($q^{(1)}_{j,d+l}$). However, this term does not contribute to the trace of $Q_{m}^{2}$ at $O(\alpha)$. Thus, we conclude that along a complete basis for the tangent space of $\Omega$ at $C_{\text{slow}}$, $D_{\Psi_{ij,kl}}\text{Tr}\left(Q_{m}^{2}\right)=0$. From this, we conclude that $\nabla_{C}\text{Tr}(Q_{m}[C_{\text{slow}}]^{2})=0$ on $\Omega$. The proof of Lemma 2 follows from application of Lemma 3. The proof of Theorem 2 then follows upon combining Lemma 1 and Lemma 2. $\square$

###### Lemma 1 .

###### Proposition 1 .

###### Lemma 2 .

###### Proposition 2 .

###### Lemma 3 .

### A.3 Supplementary Figures

Figure: Supplementary Figure 1: (Black) Average subspace angles between FCCA and PCA projections applied to Dale’s law constrained linear dynamical systems (LDS) as a function of non-normality. (Blue) Subspace angles between FCCA and PCA projections applied to firing rates derived from spiking activity driven by Dale’s Law constrained LDS. Spread around both curves indicates standard deviation taken over 20 random generations of $A$ matrices and 10 random initializations of FCCA. Panels a-d report results at projection dimension $d=4,6,8,10$, respectively, to complement the results shown in Figure 2 in the manuscript, demonstrating that non-normality drives the divergence between FCCA and PCA subspaces.

Figure: Supplementary Figure 2: Average subspace angles between $d=2$ FCCA and PCA projections applied to a switching linear dynamical system sequence as a function of non-normality. Grey spread indicates standard deviation over 20 randomly generated sequences of $A$ matrices and 10 repetitions of generated firing rates from each model sequence.

We found that PCA and FCCA identify distinct subspaces in non-normal systems. To evaluate to what degree this observation is robust to non-stationarity in the data generating process, we simulated data from a system that switched between a sequence of three $A$ matrices (still constrained to follow Dale’s law) in eq. LABEL:eq:lds_stochastic with roughly equivalent degree of non-normality. In Supplementary Figure 2, we plot FCCA/PCA subspace angles as a function of non-normality (spread taken over 20 different sequences of $A$ matrices and 10 draws of activity from each sequence). We find subspace angles to be consistently large, with only a weak dependence on non-normality. Thus, even in the case of a linear switched system [Linderman et al., ], FCCA and PCA identify distinct subspaces of dynamics.

Figure: Supplementary Figure 3: (a) Full range of average subspace angles at projection dimension $d=2$ between PCA and FCCA solutions for various T. Spread is taken over recording sessions and folds of the data within each recording session. (b) Full range of spread in average subspace angles between FCCA for $T=3$ and PCA taken across 20 initializations of FCCA and all recording sessions.

In Supplementary Figure 3, we investigate the robustness of the substantial subspace angles between FCCA and PCA observed in Figure 3a to three sources of potential variability: (i) choice of the $T$ parameter within FCCA, (ii) the dimensionality of projection, and (iii) different initializations of FCCA. In Supplementary Figure 3 a, we plot the full range of average subspace angles across recording sessions at projection dimension $d=2$ between PCA and FCCA for various choices of $T$ ($T=3$ is shown in Figure 3a). We observe that subspace angles remain consistently large ($>3\pi/8$ rads) across $T$. In Supplementary Figure 3b, we plot the full range of average subspace angles between FCCA (using $T=3$) and PCA across a range of projection dimensions. The spread in boxplots is taken across both recording sessions and twenty initializations of FCCA. We observe relatively little variability in the average subspace angles for a fixed projection dimensionality. As the projection dimension is increased, we observe the average subspace angles between FCCA and PCA decrease, from $\approx 3\pi/8$ rads to $\approx\pi/4$ rads. This is to be expected, as it is in general less likely that higher dimensional subspaces will lie completely orthogonal to each other. Overall, we conclude that FCCA and PCA subspaces are geometrically distinct in the hippocampal dataset examined.

To evaluate the robustness of FCCA’s behavioral predictions to different intializations of the algorithm, we trained linear decoders of rat position from FCCA subspaces obtained from each of twenty initializations of FCCA within each recording session. In Supplementary Figure 4, we plot the full spread in the resulting cross-validated $r^{2}$ relative to the median cross-validated $r^{2}$ as a function of projection dimension. By $d=6$, the range of spread in prediction $r^{2}$ is less than the corresponding difference between FCCA and PCA $r^{2}$. We therefore conclude that the behavioral prediction performance of FCCA is robust to the non-convexity of its objective function.

Figure: Supplementary Figure 4: Full range of variation in cross-validated position $r^{2}$ from projected FCCA activity relative to the median cross-validated $r^{2}$. Spread is taken across 20 initializations of FCCA and across all recording sessions

## Peer Reviews

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: How would Figure 3b look if compared with reduced rank regression with the relevant number of dimensions? i.e., what is the best decodability with a limited number of dimensions? To follow up on this further, are the dimensions identified by reduced rank regression also somewhat orthogonal to those identified by PCA?
soundness: 4
strengths: The paper is well-motivated and well-reasoned. The findings are very interesting and highly relevant to the neuroscience and data analysis communities. The conceptual advances are very high.
summary: This paper provides a very interesting way of thinking about relevant subspaces in recorded data. Firstly, the paper relates PCA to a specific form of feedforward control, and shows that PCA recovers the dimensions in the neural data that are most affected by an external control input. Secondly, the paper identifies a low-dimensional controller to control recorded dynamics, and defines the most controllable dimensions as the 'Feedback Controllability' components. This paper provides a conceptual advance to the field and further shows that non-normal dynamics that arise in constrained networks lead to orthogonal PCA and FCCA components, and that the feedback controllable components provide a good reconstruction of behavior.
weaknesses: The assumptions that the input is purely white noise is potentially problematic. Could the authors show, at least in simulation, that their main results hold with temporally filtered signals?

While Theorem 1 is helpful for the paper, it may not be necessary to restate it in its entirety, or might be sufficient in the Appendix.

Some key references may be missing, such as to the Henrici metric for non-normality.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: I have no problem with the math and methods developed in this study, but I have a major conceptual question about feedback control.

Conceptually, feedback control means the controller sends feedback signals back to the system (neural circuits) based on the system's outputs. Conventionally, people think the behavioral output is a (feedforward) readout of the internal state of neural circuits. Then a gap is why behavioral outputs need to send feedback to neural circuits. Further, why do behavioral outputs need even to identify the feedback controllable subspaces? Does this operation complicate the information processing or bring some actual benefits to the brain?
I can accept that the motor system and motor output (behavior) need feedback control, but I am still struggling to understand why the sensory cortex also needs that (S1 random in Table 1). I think some conceptual explanations about this are quite helpful, rather than just using this as a data analysis method to show its improved performance.
soundness: 2
strengths: 1. It seems to be novel to apply the feedback controllable subspace to analyze neural data, while the theory of feedback control (subspace) theory has been well grounded.

2. The writing of the paper is good and structure-wise. The math derivations are well laid out.
summary: This study proposes a method to identify feedback controllable subspaces from neural population responses. It shows the feedback control can identify different subspaces with feedforward control, and can explain behavioral output better.
weaknesses: 1. The data analysis seems a bit preliminary. I suggest the author could identify more intricate details of the identified subspaces from feedforward and feedback controls, and provide physical and behavioral interpretations of these subspaces if possible.

2. The equation just below Eq. 6: should the last $x$ be $x_0$?

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. The cost function for LQR (line 173) is ambiguous. If the integral is Riemann, there should be no $\frac{1}{T}$ before the integration [Burl (1998), p.283]. If it is instead an Ito integral, please indicate this and use the notation in [Jonckheere & Silverman (1983)], as this refers not to the limit but to the “limit in the mean.”

---
2. The proof of Theorem 2 (Appendix A.3) is difficult to follow due to the contradictions and typographical errors noted in the weaknesses. Once the above weaknesses are resolved, please review the proof of Theorem 2 carefully to ensure its notations are consistent with the main text. Additionally, provide references for any equations that are sourced from other works but not derived here.

---
3. [Page 7, footnote] The normal matrix can be neither symmetric nor orthogonal. $A$ is a normal matrix iff $AA’ = A’A$. One example is
$$
A=\left[\begin{array}{ccc}
1 &1 & 0 \\
0 & 1 & 1 \\
1 & 0 & 1
\end{array}\right]
\Rightarrow
AA' = A'A = \left[\begin{array}{ccc}
2 &1 & 1 \\
1 & 2 & 1 \\
1 & 1 & 2
\end{array}\right]
$$
Clearly, $A$ is neither symmetric nor orthogonal. Please fix this terminology throughout the paper.

---
4. A complete set of figures for other data types (M1 random, S1 random, and M1 maze) should be included in the appendix. Each dataset should have a figure similar to Figure 3, in addition to the information in Table 1. This is crucial to demonstrate that FCCA is applicable across multiple datasets, especially as this is the only data analysis presented in the paper.

---
**Minor problems:**
1. [Line 129] In equation (3), the dummy variable $dt$ should be after $e^{At}BB'e^{A't}$.
2. [Line 141] The volume of reachable state space is proportional to $\sqrt{\det(\Pi)}$, not $\det(\Pi)$, since the volume formula of an $n$-dimensional ellipsoid is $\frac{\pi^{n/2}}{\Gamma(n/2+1)} \prod_{k=1}^{n} r_k$ where $r_{1:n}$ are radius of the ellipsoid.
3. [Line 479] This is a typo. It should be Fig. 3b. There is no Fig. 3c.
4. [Line 325-326] “It can be shown that… to the truncated system…” Please provide references for this statement since it can be shown.
soundness: 1
strengths: - Authors develop FCCA, which computes the feedback control invariant, $\text{Tr}(PQ)$, directly through the observed neural state $x(t) \in \mathbb{R}^N$.
summary: This manuscript proposes an algorithm, Feedback Controllability Components Analysis (FCCA), to identify the low-dimensional subspace critical for feedback control. The authors claim that FCCA is a dimensionality reduction method that encodes controllability within the data. The authors apply FCCA to both simulations and data analyses.
weaknesses: 1. Throughout the paper, almost every equation lacks either a reference or a derivation, making it difficult for readers to follow. For instance:
* Section 2.1: The Lyapunov equation (3) has no reference. A reference can be [Burl (1998), p.72].
* Section 2.2: The cost functions of the Kalman filter and LQR, the dual Riccati equations (5) and (6), and their associated parameter equations for $Q$ and $P$ lack references. References can be, for example, [Burl (1998), p.243 & p.283] and [Jonckheere & Silverman (1983)].
* Section 2.3: The final FCCA equation (11) has no derivation details. Considering it's the main result of this paper, the derivation should be step-by-step. Additionally, the authors do not provide pseudocode for FCCA, making the algorithm challenging to follow.

---
2. The notation throughout the paper is inconsistent, with numerous typos. For example, in Section 2.2, the term $QC’CQ$ in equation (5) should appear only when the observation model is $y(t) = Cx(t) + v(t)$, where $v(t)$ is Gaussian white noise with covariance $I_d$. Sources such as [Burl (1998), p.243] and [Jonckheere & Silverman (1983)] include $v(t)$, leading their Riccati equations to include $QC’CQ$. However, [Ljung & Kailath (1976)] exclude $v(t)$, and thus their Riccati equation does not include $QC’CQ$. This creates a contradiction:
* The $QC’CQ$ term seems necessary in deriving FCCA in equation (9) because it is dual to the LQR Riccati equation (6).
* However, since the neural state $x(t)$ represents the true observation and $y(t) = Cx(t)$ is only a low-dimensional projection, there should be no noise term $v(t)$.
It is unclear how the authors intend to resolve this contradiction. This further underscores the importance of providing step-by-step derivations from the base model to ensure clarity and avoid such contradictions, making it easier for readers to follow. The authors should clarify their assumptions about the observation model and provide a step-by-step derivation showing how they arrive at equation (5) from their model equation (2), particularly explaining the presence or absence of the $QC'CQ$ term.
---
3. [Kashima (2016)] has shown that the feedforward cost matrix (i.e., the controllability Gramian matrix $\Pi$) is equivalent to the data covariance. Additionally, Section 2.2 on feedback controllability closely follows the results in [Jonckheere & Silverman (1983)], including Theorem 1 and the similarity invariant matrix $QP$. Therefore, the only novel contribution of this paper is the derivation of FCCA in Section 2.3. However, this derivation contains serious typographical errors and lacks sufficient detail.

The first issue with FCCA is that it should be $x_b(t) = \Pi x_a(t)$ rather than $x(t) = \Pi x_a(t)$ in line 239, which would allow the authors to derive the equation in line 241 from equation (8). Although this may seem like a minor typo, it is critical because $x_a(t)$ is no longer connected to $x(t)$. Since $x(t)$ represents the observed neural state rather than a latent state, it is unclear how the covariance and cross-covariance of $x_a(t)$ could be estimated from $x(t)$ as required in equation (11), the FCCA formula. The authors should provide a detailed explanation of how they connect $x_a(t)$ to the observed neural state $x(t)$, and how this connection allows for the estimation of covariances $\tilde{P}$ in equation (11).

---
4. The second issue with FCCA is that the Riccati equation (9) for $\tilde{P}$ is not the same as the Riccati equation (6) for $P$. It is unclear why the authors equate $\text{Tr}(QP)$ with $\text{Tr}(Q\tilde{P})$. A detailed derivation linking these two expressions is necessary. Additionally, it is not explained why equation (9) aligns with the Riccati equation associated with the modified LQR cost function (10). Since $C’C$ is replaced by $\Pi^{-1}BB’\Pi^{-1}$ in (10), the final LQR Riccati equation should exclude $C’C$, which contradicts equation (9). A detailed derivation of this step is also required. The authors should provide a step-by-step derivation showing how equation (9) can be transformed into equation (6), so they're equivalent, and how it relates to the modified LQR cost function in equation (10).

---
5. The third issue with FCCA is that equation (11) lacks a derivation. This should be derived step-by-step, as it represents the main result of the paper. Furthermore, some details in equation (11) appear questionable. For instance, following my previous point that $\Pi x_a(t) = x_b(t)$ rather than $x(t)$, the inverse matrix $\Sigma_T^{-1}(C)$ of $\tilde{P}$ in (11) should be based on $x_b(t)$, not $x(t)$. Therefore, this formula appears to be incorrect. These issues need to be clarified through a detailed derivation. The authors should provide a complete, step-by-step derivation of equation (11), starting from the basic assumptions and clearly stating any approximations or simplifications made along the way.

---
6. Beyond the theoretical issues mentioned, a more fundamental question arises: Why not simply fit the parameters $A$ and $B$ of the linear dynamical model in equation (2) directly and solve the dual Riccati equations (5) and (6) iteratively? This would yield the Gramian matrix $\Pi$, the Kalman error covariance $Q$, and the LQR cost matrix $P$. This approach is straightforward since $x(t)$ is the observed neural signal, not a hidden latent state, making it feasible to fit $A$ and $B$ through simple linear regression. Is there any practical advantage to using FCCA equation (11) compared to this more direct approach? A detailed comparison between the FCCA method and the more direct approach, including computational complexity, accuracy, and any other relevant factors will be helpful.

---
7. The simulations do not validate the correctness of FCCA. To achieve this, the authors should simulate a linear dynamical system (which appears to be the LDS model in Figure 2), compute $Q$ and $P$ using the dual Riccati equations (5) and (6), and find $C$ by minimizing $\text{Tr}(QP)$. Then, they should compute $Q$ and $\tilde{P}$ using equation (11) and confirm that these values match $Q$ and $P$ from equations (5) and (6). Finally, the optimal $C$ derived from (11) must align with the original $C$ that minimizes $\text{Tr}(QP)$. This verification should be straightforward since all system parameters are precisely defined in the simulations. Please perform this validation and include the results in a new figure or table in the paper.

---
8. The simulations also fail to validate the correctness of Theorem 2. The authors should simulate the case where $B = I_N$ and $A = A’$, which would make the angle between the PCA and FCCA subspaces equal to zero. In other words, Figure 2a should include the case where $\|AA'-A'A\|_F=0$.


---
Overall, while the overall idea of finding the subspace with feedback controllability is nice, the derivation of FCCA is incomplete and unconvincing and novelty is unclear compared with prior work. Furthermore, the simulations do not properly validate the method's correctness and the data analyses are limited and unconvincing for demonstrating that the method works.


**References**
Burl, J. B. (1998). Linear optimal control: H (2) and H (Infinity) methods. Addison-Wesley Longman Publishing Co., Inc..

---

confidence: 2
contribution: 3
initial_rating: 3
presentation: 1
questions: - Is the controllability Gramian $\Pi$ equal to the covariance of the stationary distribution of $x(t)$ defined in eq. (2)? I think this is a more intuitive description that could be stated earlier.
- Notation in section 2.2: on line 182, $Q$ is set equal to the (arg)-min over probability distributions. It appears this means that the optimal distribution is Gaussian with covariance matrix $Q$. Are these minimizations equivalent since the optimal distribution is a mean-zero Gaussian? It might help to state this.
- Line 188: What is the meaning of "$P$ encodes the regulation cost incurred for varying initial conditions"?
- Line 222: What is causal vs acausal Kalman filter? Does acausal mean the Kalman filter applied to the backward process $x_b(t)$?
- Theorem 2: If $A$ is symmetric, then isn't $Re(\lambda(A))=\lambda(A)$? If so, I'd suggest stating this.
soundness: 3
strengths: Identifying relevant low-dimensional projections of high-dimensional neural data is a highly relevant topic. This is evidenced by the numerous methods that have been proposed in the literature. This paper proposes a compelling methodology that is based on the controllability of the neural state space. They build on a well-developed theory from the control theory community that is not well-known in the neuroscience community.
summary: This paper proposes a dimensionality reduction method called **Feedback Controllability Components Analysis (FCCA)**. The method identifies the subspace of a linear dynamical system that is maximizes a measure of feedback controllability, which is defined in terms of the optimal Kalman filter. Specifically, the *feedback controllability* is the trace of the product of two positive definite matrices: one matrix encoding the covariance of the estimation error and another matrix encoding how the optimal control cost depends on initial condition. They compare FCCA with PCA and show that the angle between the subspaces depends on the non-normality of the dynamics matrix. They also apply their method to neural recordings from rat hippocampus, macaque primary motor cortex and primary somatosensory cortex. In each case, they find that the FCCA projection is more predictive of the animal's behavior than the PCA projection.
weaknesses: - I found the paper to be dense and required significant effort to understand (at least for me, the main reason for my presentation score of 1 and my overall recommendation for rejection). One of the stated goals of this paper is to highlight important ideas from the control theory community that are relevant to the neuroscience community. I do not think this paper lays out the ideas with sufficient clarity for the theoretical neuroscience community (my confidence score is 2 because I don't think I fully understand the paper). While some of this is perhaps unavoidable due to the material, I think the presentation could be substantially improved with an intuitive figure explaining the concepts in this paper. It would be very helpful to have an illustrative example (e.g. in 2D) of the relationship between $A$, $B$ and the Gramian $\Pi$ and how PCA and feedback controllability subspaces differ. A large portion of the computational neuroscience community is quite familiar with 2D linear dynamical systems and I think this paper misses an opportunity to connect the results to prior understanding in the community.
- FCCA is only compared with PCA in section 4. There are a multitude of methods for extracting subspaces beyond PCA; e.g., slow feature analysis, GFPA, LFADS, etc. Such comparisons seem important if the goal is to encourage practitioners to use FCCA.