Closed-form control with spike coding networks

We mathematically formalize the link between optimal control of dynamical systems and the SCN framework. Particularly, we analytically derive the spiking equivalent of the linear–quadratic–Gaussian (LQG) control problem – the combination of a linear quadratic estimator (Kalman filter) and a linear–quadratic regulator (LQR) controller (Fig. 2B and C). We show that our proposed networks ($i$) accurately estimate and control well-known dynamical systems, achieving similar performance as their non-spiking counterpart, and, importantly, ($ii$) they preserve the irregular and sparse spiking patterns and robustness to neural silencing.

Developing efficient spiking neural controllers with fully derivable and interpretable connectivity is especially relevant for industrial applications [6], robotics [9] and machine intelligence in general [12].

Here, we non-exhaustively review the main lines of research that followed a similar approach as us. For recent reviews on SNNs related to robotics and control see [8, 22]. Our work is based on [15], which originally showed how any linear dynamical system can be analytically implemented in an SCN. While there have been several follow-up papers on extending the framework for general computations [19, 23, 18], control has been less well-studied. In [21], spike coding networks were extended by adding learning rules for connectivity weights. This allowed them to perform forward prediction of the dynamical system state. As an additional result, they showed how to control a pendulum by using the network to simulate multiple future trajectories of the pendulum with different control policies. However, the control algorithm is not encoded within the network and does not provide robust estimation. In [24, 25], the authors proposed an analytical SCN-inspired framework for control. They mainly focused on producing the correct control signal, and derived network connectivities accordingly. The framework requires fully observable systems, and they were not able to provide a general, clear and simple mathematical model that is in line with both control and neuroscience communities. Furthermore, they did not investigate the beneficial properties of the spiking patterns generated by the coordination of the neurons. Here, we will provide a simpler analytical derivation more in line with existing SCN theory, derive networks for both estimation and control of partially observable systems, and investigate the robustness properties of the control and networks in more detail. Overall, our approach takes inspiration from the broader SCN literature to provide a unified mathematical framework to analytically compute optimal control in SNNs.

We have a system with state $\mathbf{x}(t)$ we would like to control with a spiking signal ($\mathbf{s}(t)$) emitted from a recurrent SNN (Fig. 1). The SNN is provided with some incomplete measurements of the system state $\mathbf{y}(t)$ and a target state $\mathbf{z}(t)$, and has to generate a control signal $\mathbf{u}(t)$.

The spiking patterns are generated according to some underlying voltage dynamics. There are many such models of varying complexities, but for most practical applications networks of leaky integrate-and-fire (LIF) neurons suffice. A network of $N$ such neurons is then defined by

where $\mathbf{v}(t)\in\mathbb{R}^{N}$ is the vector of neural voltages, $\mathbf{F}\in\mathbb{R}^{K\times N}$ are some input weights, $\mathbf{c}(t)\in\mathbb{R}^{K}$ is some $K$-dimensional input, $\mathbf{s}(t)\in\mathbb{R}^{N}$ are the emitted spikes, $\mathbf{r}(t)\in\mathbb{R}^{N}$ are the filtered spike-trains¹¹1such that $\dot{\mathbf{r}}=-\lambda\mathbf{r}+\mathbf{s}$, see Fig. 1., and $\mathbf{\eta}_{v}\in\mathbb{R}^{N}$ corresponds to independent voltage noise. Finally, there are fast synapses which affect the post-synaptic voltage instantaneously following a spike (through fast recurrent connections $\mathbf{\Omega}_{f}\in\mathbb{R}^{N\times N}$), and slow synapses which cause an initially slow change in voltage (through the slow recurrent connections $\mathbf{\Omega}_{s}\in\mathbb{R}^{N\times N}$).

Whenever a given neuron’s voltage $v_{i}$ reaches a threshold $T_{i}$, that neuron will emit a spike at that time (resulting in spike-trains $s_{i}(t)=\sum_{j}\delta(t-t_{j})$). The voltage will then be reset to some resting value (here implemented through the diagonal elements of $\mathbf{\Omega}$).

How should the recurrent and input weights of the network be set, such that the output signal is as optimal as possible? One solution is to train the network using a cost function that enforces, for instance, sparsity [14]. Here we take another route by first defining the solution according to classical control theory, and then translating this solution into the spiking network parameters by using the neuroscience theory of Spike Coding Networks (SCN) [15, 26].

For the sake of completeness, in this section we will give a brief overview of the original SCN framework (as originally proposed in [15]) for implementing linear dynamical systems in an SNN — for an extensive mathematical derivation of the SCN we refer the reader to [17, 18].

Consider a $K$-dimensional linear dynamical system in state-space representation (Fig. 2A.i):

where $\mathbf{x}(t)\in\mathbb{R}^{K}$ is the state vector, $\mathbf{A}\in\mathbb{R}^{K\times K}$ is the system matrix conveying the dynamics of the system, $\mathbf{u}(t)\in\mathbb{R}^{P}$ is the input vector through which we can control the system, and $\mathbf{B}\in\mathbb{R}^{K\times P}$ is the input matrix. $\boldsymbol{\eta}_{d}\in\mathbb{R}^{K}$ reflects internal disturbances, given by a zero-mean Gaussian process with co-variance $\mathbf{\Sigma}_{d}$. $\mathbf{y}(t)\in\mathbb{R}^{Q}$ are observations about $\mathbf{x}(t)$, where $\mathbf{C}\in\mathbb{R}^{Q\times K}$ is the observability matrix. $\boldsymbol{\eta}_{n}$ is sensor noise, given by a zero-mean Gaussian process with co-variance $\mathbf{\Sigma}_{n}$.

We will generally assume as known the measurement vector $\mathbf{y}$, the system matrix $\mathbf{A}$, input matrix $\mathbf{B}$, the input vector $\mathbf{u}$, the observability matrix $\mathbf{C}$, and co-variances $\mathbf{\Sigma}_{d}$ and $\mathbf{\Sigma}_{n}$.

Our goal is to control $\mathbf{x}$ to some reference state $\mathbf{z}\in\mathbb{R}^{K}$. To do this we must be able to estimate the state $\mathbf{x}$ from the observations $\mathbf{y}$ (estimation problem), and find the best control signal $\mathbf{u}$ to do so (optimal control problem). In control theory this is largely solved, and we will now consider how to combine the resulting solutions with SCN theory to generate the correct u as an output of an SNN.

Full-state estimates $\mathbf{\hat{x}}$ given a noisy and incomplete measurements in linear systems can be provided by a Kalman filter [27], which optimally balances an internal (predicted) state estimate to a noisy and/or partially observable external measurement (Fig. 2B.i). The Kalman filter in continuous time is given by a dynamical system of the form

The Kalman filter gain matrix, $\mathbf{K}_{f}\in\mathbb{R}^{K\times Q}$, is applied to an error between the observations of the dynamical system and the Kalman filter’s own internal estimate, $\mathbf{\hat{y}=C\hat{x}}$. For fully identified systems $\mathbf{K}_{f}$ can be found by solving an algebraic Riccati equation [27]. Assuming that $\mathbf{K}_{f}$ is known, we can directly implement the optimal filtering under the SCN framework (by following Section II-B2) resulting in the following voltage update rule:

where, on top of the previously introduced slow and fast connectivity, we now also have input (control) connections $\mathbf{F}_{i}=\mathbf{D}^{\top}\mathbf{B}$ mapping the control signal $\mathbf{u}$ to the neurons, recurrent and feed-forward “Kalman filter” connections ($\mathbf{\Omega}_{k}=-\mathbf{D}^{\top}\mathbf{K}_{f}\mathbf{CD}$ and $\mathbf{F}_{k}=\mathbf{D}^{\top}\mathbf{K}_{f}$). The recurrent connections essentially read-out the internal state of the system, apply both the observability matrix and the Kalman filter gain matrix to predict its evolution, and maps this back into the network. The feed-forward connections take the partially observable state, $\mathbf{y}$, and applies the Kalman filter gain matrix. A more detailed derivation can be found in the appendix (section VIII). Now, the entire SCN, including all of its connections, represents a Kalman filter (Fig. 2B.ii). The network internally simulates a linear dynamical system of the form $\mathbf{\dot{\hat{x}}=A\hat{x}+Bu}$, but also corrects its own estimate according to the external input $\mathbf{y}$.

Now that we have optimal state-estimation we can extend the framework to optimal control by deriving the spiking version of a Linear-Quadratic Regulator (LQR) controller. Given a known system of the form in Eq. (5), this controller produces the optimal $\mathbf{u}$ signal to minimize the squared error between the state $\mathbf{x}$ and a target $\mathbf{z}$. This gives rise to a linear control law of the form

where $\mathbf{K_{c}}\in\mathbb{R}^{P\times K}$ is the LQR gain matrix, which can be found by solving an algebraic Riccati equation, given assumptions on the cost of state deviations and actuation [27].

In order to control a partially observable and/or noisy dynamical system, we can combine an LQR controller with a Kalman filter (see Fig. 2C.i). The Kalman filter should then be aware of the specific form used for the control, and we can swap out $\mathbf{u}$ for the LQR control law (Eq. (9) in Eq. (7)), resulting in a dynamical system of the form

We will refer throughout to the combination of a Kalman filter and LQR control as the ideal controller. We can now use the same method as in the previous section to implement this in a spiking network.

To be able to read out the full control signal, we additionally encode the reference $\mathbf{z}$ into the network with a new set of fast connections (similar to Eq. (2)), which allows us to read out $\mathbf{u}$ from the internal estimates $\mathbf{\hat{x}}$ and $\mathbf{\hat{z}}$ of the network. For this we define two sets of decoding weights: $\mathbf{D}_{\mathbf{x}}$ for $\mathbf{\hat{x}}$, and $\mathbf{D}_{\mathbf{z}}$ for $\mathbf{\hat{z}}$.

Extending the SCN defined in Eq. (8) with LQR control and an internal encoding of $\mathbf{z}$, we get the final voltage update rule

where we have indicated the parts of the connectivity that are tracking the external system, the effect of the control on the system, the Kalman filter updates, and the representation of the target signal. For the internal estimate $\mathbf{\hat{x}}$, we now have two sets of recurrent “control” and “target” connections, which represent the LQR controller within the Kalman filter SCN (see Fig. 2C.ii). The recurrent control connectivity, $\mathbf{\Omega}_{c}=-\mathbf{D}_{\mathbf{x}}^{\top}\mathbf{B}\mathbf{K}_{c}\mathbf{D}_{\mathbf{x}}$, decodes $\mathbf{\hat{x}}$ from the internal state of the SCN using $\mathbf{D}_{\mathbf{x}}$, and applies the LQR gain matrix $\mathbf{K}_{c}$. The result is transformed using $\mathbf{B}$, and encoded back into the network. The recurrent target connectivity, $\mathbf{\Omega}_{z}=\mathbf{D}_{\mathbf{x}}^{\top}\mathbf{B}\mathbf{K}_{c}\mathbf{D}_{\mathbf{z}}$ does something similar, but first reads out $\mathbf{\hat{z}}$ using $\mathbf{D}_{\mathbf{z}}$. The reference signal, $\mathbf{z}$, is encoded into the network using the same method as defined in Eq. (8). A more detailed derivation for all these connections can be found in the appendix (section VIII). Note that while there now appear to be a large set of separate recurrent connections, all the slow connectivities ($\mathbf{\Omega}_{s}$, $\mathbf{\Omega}_{c}$, $\mathbf{\Omega}_{z}$, $\mathbf{\Omega}_{k}$) can be grouped together in a single connection matrix, as well as both the fast connectivities ($\mathbf{\Omega}_{f}$, $\mathbf{\Omega}^{\mathbf{z}}_{f}$).

We now have an SCN which implements the entirety of $\mathbf{\dot{\hat{x}}=A\hat{x}+B}\mathbf{K}_{c}(\mathbf{\hat{x}}-\mathbf{\hat{z}})+\mathbf{K}_{f}\mathbf{(\hat{y}-y)}$, but also internally keeps an estimate of the desired state or reference $\mathbf{z}$. We can obtain the output of the internal LQR controller in the SCN by defining the following new set of decoding weights: $\mathbf{D_{u}}=-\mathbf{K}_{c}(\mathbf{D}_{\mathbf{x}}-\mathbf{D}_{\mathbf{z}})$. The control signal can then directly be read-out from the neural activities as $\mathbf{u}=\mathbf{D_{u}}\mathbf{r}$. This $\mathbf{u}$ can be applied to an external dynamical system to control it, all whilst the SCN keeps an internal estimate of the external dynamical system. Hence, we have finally arrived at a fully derived recurrent SNN which allows us to track, replicate and control an external dynamical system.

We compared our proposed SCN estimation and control with their non spiking counterparts for the spring-mass-damper system (Fig. 3A). The baselines comparison are a Kalman filter for estimation and LQG for control, which we refer to as the ”Idealized” estimator and controller. Furthermore, we evaluated the system robustness against input-noise, neural silencing and external perturbations. The SMD system has two dynamical variables: position $x$ ($=x_{1}$) and velocity $v$ ($=x_{2}=\dot{x}_{1}$). Our instance of the SMD system has an internal disturbance $\mathbf{\eta}_{d}$, and it outputs a partially observable state which measures only the position $x$ with sensor noise $\mathbf{\eta}_{n}$.

We evaluated our proposed SCN controller (i.e., linear-quadratic-Gaussian spiking control) on the nonlinear cartpole system (see Fig. 5A). The cartpole has four dynamical variables: the position of the cart, $x$ ($=x_{1}$), the velocity of the cart $v$ ($=x_{2}=\dot{x}_{1}$), the position of the pole, $\theta$ ($=x_{3}$), and the angular velocity of the pole, $\omega$ ($=x_{4}=\dot{x}_{3}$). Just like the SMD, our instance of the cartpole has an internal disturbance $\mathbf{\eta}_{d}$, and it outputs a partially observable state consisting of only the position of the cart, $x$, with sensor noise $\mathbf{\eta}_{n}$. Note that the cartpole system is a nonlinear dynamical system, but both our ideal and spiking controllers assume linear dynamics. To produce a meaningful control signal, we use linearized dynamics for the internal estimate. In our case, we linearized around the up-position of the pole ($\theta=\pi$).

Fig. 5B compares a cartpole system controlled through an idealized controller (red) to a cartpole controlled through an SCN controller (green). On top of that, we show the reference signal (blue, dashed), which is a stair-wise increase of the cart position, $x$. The task of the controllers is then to move the cart whilst keeping the pole upright. In the figure, we show the position ($x$) of the cart and the angle of the pole ($\theta$). We see that both controllers follow the reference signal almost perfectly, and most importantly, keep the pole upright.

The spike plots show that especially with a large amount of neurons, the spiking is irregular and sparse. Increased spiking activity is observed when the reference signal demands a change of the state of the dynamical system.

For each network, the decoding weights ($\mathbf{D}$) are sampled from a normal distribution, with each column normalized to have norm 0.1, except for norm 0.01 in Fig. 5B and norm 1 in Fig. 6. 20 and 50 neurons were used in Fig. 3B and C respectively, 50 neurons were used in Fig. 4 and Fig. 6, and 100 neurons in Fig. 5B. The SMD parameters were $m=3$, $k=5$ and $c=0.5$ for Fig. 3B and Fig. 4, while $m=20$, $k=6$ and $c=2$ were used for Fig. 3C and Fig. 6. Fig. 3B and Fig. 6 used $\mathbf{\Sigma}_{d}=\mathbf{\Sigma}_{n}=0.001$. Fig. 3C used $\mathbf{\Sigma}_{d}=\mathbf{\Sigma}_{n}=0.1$. Fig. 4 used $\mathbf{\Sigma}_{d}=0.001$. The cartpole system parameters were $m=1$, $M=5$, $L=2$, $g=-10$, and $d=1$, with $\mathbf{\Sigma}_{d}=\mathbf{\Sigma}_{n}=1\text{e-}7$. All networks used $\eta_{V}=1\text{e-}5$ and $\lambda=0.1$, except for Fig. 6, which used $\eta_{V}=1\text{e-}6$ and varying $\lambda$. Forward Euler was used throughout, with a timestep of 1e-3s, except for 1e-4s in Fig. 4, Fig. 5B and Fig. 6. For control, $R=1\text{e-}2$. For control of the SMD system, the matrix $\mathbf{Q}\in\mathbb{R}^{K\times K}$ prioritizes $x_{1}$ with weight 10, and $x_{2}$ with weight 1. For cartpole control, $x_{3}$ is weighted with 10 and the rest with 1.

For the continuous definition of the Kalman filter (Eq. 7), the above results in voltage dynamics given by

We can then replace $\hat{\mathbf{x}}=\mathbf{D}\mathbf{r}$, and group the different terms together to get the final voltage dynamics to implement a Kalman filter in a recurrent SNN given by Eq. 8.

For the continuous definition of a Kalman filter including LQR control we have equation 10. We can do the same as above and directly implement this in a recurrent SNN, but we also additionally need to represent the target state $\mathbf{z}$ in the network (as explained in the main text). This together results in the voltage update rule

We can then replace $\hat{\mathbf{x}}=\mathbf{D}\mathbf{r}$, and group the different terms together the get the final voltage dynamics to implement a Kalman filter plus LQR controller in a recurrent SNN given by Eq. 11.