# Learning to Time-Decode in Spiking Neural Networks Through the Information
Bottleneck

Nicolas Skatchkovsky
KCLIP lab, Dept. of Eng.
King’s College London
London, UK
<EMAIL_ADDRESS>
&Osvaldo Simeone
KCLIP lab, Dept. of Eng.
King’s College London
London, UK
<EMAIL_ADDRESS>
&Hyeryung Jang
ION group, Dept. of A.I.
Dongguk University
Seoul, Korea
<EMAIL_ADDRESS>

###### Abstract

One of the key challenges in training Spiking Neural Networks (SNNs) is that
target outputs typically come in the form of natural signals, such as labels
for classification or images for generative models, and need to be encoded
into spikes. This is done by handcrafting target spiking signals, which in
turn implicitly fixes the mechanisms used to decode spikes into natural
signals, e.g., rate decoding. The arbitrary choice of target signals and
decoding rule generally impairs the capacity of the SNN to encode and process
information in the timing of spikes. To address this problem, this work
introduces a hybrid variational autoencoder architecture, consisting of an
encoding SNN and a decoding Artificial Neural Network (ANN). The role of the
decoding ANN is to learn how to best convert the spiking signals output by the
SNN into the target natural signal. A novel end-to-end learning rule is
introduced that optimizes a directed information bottleneck training criterion
via surrogate gradients. We demonstrate the applicability of the technique in
an experimental settings on various tasks, including real-life datasets.

##  1 Introduction

While traditional Artificial Neural Networks (ANNs) implement static rate-
based neurons, Spiking Neural Networks (SNNs) incorporate dynamic spike-based
neuronal models that are closer to their biological counterparts. In different
communities, SNNs are used either as models for the operation of biological
brains or as efficient solutions for cognitive tasks. The efficiency of SNNs,
when implemented on dedicated hardware, hinges on the high capacity, low
latency, and high noise tolerance of information encoding in the timing of
spikes [1]. Recent results leveraging novel prototype chips confirm the
potential order-of-magnitude gains in terms of latency- and energy-to-accuracy
metrics for several optimization and inference problems [2]. In contrast, the
landscape of training algorithms for SNNs is still quite fractured [2], and a
key open question is how to derive learning rules that can optimally leverage
the time encoding capabilities of spiking neurons. This is the focus of this
work.

Figure 1: (a) Illustration of an SNN trained for image classification with
first-to-spike decoding. (b) Illustration of an SNN trained for image
classification with rate decoding. (c) Illustration of an SG scheme [3]: At
each time-step, updates for neuron $i$ are computed using local eligibility
traces ${\bm{e}}_{i,t}$, and per-neuron learning signals $L_{i,t}$, obtained
from errors computed based on handcrafted target signals (gray). (d)
Illustration of ML learning for SRM with stochastic threshold [4]: Updates for
readout neurons (green) are obtained with local eligibility traces based on
handcrafted target spiking signal (gray), while updates for hidden neurons
(red) make use of a common learning signal $\ell_{t}$ based on the same target
signal. (e) Illustration of the proposed hybrid variational autoencoder: An
encoding SNN processes inputs from a spiking image source, such as a DVS
camera, to produce a spiking latent representation. A decoding ANN learns to
decode this representation with the aim of reconstructing a natural target
signal, such as an image, corresponding to the spiking input of the SNN. SNN
and ANN are jointly trained under a directed information bottleneck criterion.
Accordingly, training of the ANN is based on backpropagation, while training
of the SNN uses per-neuron learning signals and common feedback from both
networks.

#### Learning how to decode.

Consider the problem of training an SNN – a network of spiking neurons. To fix
the ideas, as illustrated in Fig. 1(a)-(b), say that we wish to train an SNN
to classify inputs from a neuromorphic sensor, such as a DVS camera [5], which
are in the form of a collection of spiking signals. We hence have a spiking
signal as input and a natural signal (the label) as output. As another example
in Fig. 1(e), we may wish to “naturalize” the output of a DVS camera to
produce the corresponding natural image [6]. The direct training of an SNN
hinges on the definition of a loss function operating on the output of the
SNN. This is typically done by converting the target natural signals into
spiking signals that serve as the desired output of the SNN. As seen in Fig.
1(a)-(b), one could for instance assign a large-rate spiking signal or an
earlier spike to the readout neuron representing the correct class. The loss
function then measures the discrepancy between actual output spiking signals
and the corresponding spiking targets.

The approach outlined in the previous paragraph has an important problem: By
selecting handcrafted target spiking signals, one is effectively choosing a
priori the type of information encoding into spikes that the SNN is allowed to
produce. This constrains the capacity of the SNN to fully leverage the time
encoding capabilities of spiking signals. For instance, the target signals
outlined in Fig. 1(a) correspond to assuming first-to-spike decoding, while
those in Fig. 1(b) to rate decoding [7]. This paper proposes a method that
allows the SNN to automatically learn how to best encode the input into
spiking signals, without constraining a priori the decoding strategy.

To this end, we introduce the hybrid SNN-ANN autoencoder architecture shown in
Fig. 1(e) that consists of the encoding SNN to be trained alongside a decoding
ANN. The role of the decoding ANN is to learn how to best convert the spiking
signals output by the SNN into the target natural signal. Encoding SNN and
decoding ANN are trained end-to-end so as to ensure that the SNN produces
informative spike-based representations that can be efficiently decoded into
the desired natural signals.

#### Directed information bottleneck.

How to jointly train the hybrid SNN-ANN autoencoder? On the one hand, we wish
the spiking output of the SNN to be informative about the target natural
signal – which we will refer to as _reference signal_ – and, on the other, we
would like the spiking output of the SNN to remove all extraneous information
present in the spiking input that is not informative about the reference
signal. To capture this dual requirement, we adopt an _information bottleneck_
(IB) formulation [8], and view the ANN as the inference network in a
_variational autoencoder_ system [9].

In applying the _variational IB_ methodology to the hybrid SNN-ANN autoencoder
in Fig. 1(e), we make two key contributions. 1) Given the inherently dynamic
operation of the SNN, we generalize the IB criterion to apply to directed
information metrics [10], in lieu of conventional mutual information. We refer
to the resulting approach as variational directed IB (VDIB). 2) In order to
enable the optimization of information-theoretic metrics, such as the IB, we
introduce an SNN model that encompasses _deterministic hidden neurons and
probabilistic readout neurons_. In so doing, we combine two strands of
research on training SNNs: The first line focuses purely on deterministic
models, viewing an SNN as a Recurrent Neural Network (RNN) and deriving
approximations of backpropagation-through-time (BPTT) through surrogate
gradient (SG) techniques [11, 12, 13, 14, 15, 3]; while the second assumes
probabilistic neurons only, deriving approximations of maximum likelihood (ML)
learning. Additional discussion on related work can be found in the next
section.

##  2 Related works

#### Training SNNs.

In recent years, SG techniques for training SNNs have gained significant
momentum, becoming the state-of-the-art and de-facto eclipsing other
techniques. As reviewed in [12], these techniques rely on surrogate
derivatives to replace the ill-defined derivative of the step activation
function of spiking neurons. The SGs are combined with truncated and
approximated versions of BPTT, including the use of per-layer local error
signals [13] and per-neuron learning signals defined using direct feedback
alignment [15, 3]. Learning rules for probabilistic SNNs are typically derived
by tackling the ML problem, and result in updates that rely on random sampling
and common – rather than per-neuron – learning signals. The approach was shown
to generalize to networks of winner-take-all neurons [16]. As illustrated in
Fig. 1, the proposed approach builds on both models to enable the use of SG
methods for the optimization of probabilistic learning criterion.

#### Hybrid SNN-ANN architectures.

A hybrid SNN-ANN architecture has been recently proposed in [17] for gesture
similarity analysis, based on the variational autoencoder (VAE) framework.
Although yielding a similar objective, the encoder in [17] is not fully
spiking since the ANN receives the membrane potentials of the spiking neurons
in the readout layer of the SNN. The model proposed in our work circumvents
this limitation.

#### Information bottleneck.

The IB technique originated two decades ago [8] in the context of a rate-
distortion problem formation assuming discrete variables. The approach has
regained interest in recent years, first as a means to explain the performance
of ANNs [18], and then as a training criterion in [9]. We extend the IB by
considering temporal signals under causality constraints, hereby replacing the
mutual information with directed information [10]. The IB principle has been
previously used in the SNN literature in the case of a single neuron [19, 20].
Finally, it has been used in theoretical neuroscience to explain processing in
sensory neurons [21] and as a unifying principle for the predictive and
efficient coding theories. Our model encompasses the framework introduced in
[21].

##  3 Problem definition

#### Notations.

Consider three jointly distributed random vectors ${\bm{x}},{\bm{y}}$, and
${\bm{z}}$ with distribution $p({\bm{x}},{\bm{y}},{\bm{z}})$. The conditional
mutual information (MI) between ${\bm{x}}$ and ${\bm{y}}$ given ${\bm{z}}$ is
defined as
$I({\bm{x}};{\bm{y}}|{\bm{z}})=\mathrm{E}_{p({\bm{x}},{\bm{y}},{\bm{z}})}\Big{[}\log\Big{(}\frac{p({\bm{x}},{\bm{y}}|{\bm{z}})}{p({\bm{x}}|{\bm{z}})p({\bm{y}}|{\bm{z}})}\Big{)}\Big{]}$,
where $\mathrm{E}_{p}[\cdot]$ represents the expectation with respect to the
distribution $p$. The Kullback-Leibler (KL) divergence between two probability
distributions $p({\bm{x}})$ and $q({\bm{x}})$ defined on the same probability
space is
$\text{KL}(p({\bm{x}})||q({\bm{x}}))=\mathrm{E}_{p({\bm{x}})}\Big{[}\log\Big{(}\frac{p({\bm{x}})}{q({\bm{x}})}\Big{)}\Big{]}$.

For any two jointly distributed random vectors
${\bm{a}}=({\bm{a}}_{1},\dots,{\bm{a}}_{T})$ and
${\bm{b}}=({\bm{b}}_{1},\dots,{\bm{b}}_{T})$ from time $1$ to $T$, the
$\tau$-causally conditioned distribution of ${\bm{a}}$ given ${\bm{b}}$ is
defined as

$\displaystyle
p^{(\tau)}({\bm{a}}||{\bm{b}})=\prod_{t=1}^{T}p({\bm{a}}_{t}|{\bm{a}}^{t-1},{\bm{b}}_{t-\tau}^{t}),$
(1)

where we denote as
${\bm{a}}_{t_{1}}^{t_{2}}=({\bm{a}}_{t_{1}+1},\dots,{\bm{a}}_{t_{2}})$ the
overall random vector from time $t_{1}+1$ to $t_{2}$, and write
$a^{t}=a_{1}^{t}$; This distribution (1) captures the causal dependence of
sequence ${\bm{a}}$ on the last $\tau$ samples of ${\bm{b}}$. If $\tau=T$, we
recover the causally conditioned distribution [10].

We also define the $\tau$-directed information with $\tau\geq 0$ as

$\displaystyle
I^{(\tau)}({\bm{a}}\rightarrow{\bm{b}})=\sum_{t=1}^{T}I({\bm{a}}^{t}_{t-\tau};{\bm{b}}_{t}|{\bm{b}}^{t-1}).$
(2)

This metric quantifies the causal statistical dependence of the last $\tau$
samples from sequence ${\bm{a}}$ to ${\bm{b}}$. With $\tau=T$, this recovers
the standard directed information [10]. We denote as $\ast$ the convolution
operator $f_{t}\ast g_{t}=\sum_{\delta\geq 0}f_{\delta}g_{t-\delta}$.

#### Problem definition.

We study the hybrid SNN-ANN autoencoder architecture in Fig. 1 with the aim of
training an SNN encoder to produce a time-encoded spiking representation
${\bm{y}}$ of the spiking input ${\bm{x}}$ that is informative about a natural
reference signal ${\bm{r}}$. Throughout, a spiking signal is a binary vector
sequence, with “$1$” representing a spike, and a natural signal is a real-
valued vector sequence. In the prototypical application in Fig. 1(e), the
input ${\bm{x}}$ is obtained from a neuromorphic sensor, such as a DVS camera
[5], and the reference signal ${\bm{r}}$ represents the corresponding target
natural signal, such as a label or the original natural image. The spiking
representation ${\bm{y}}$ produced by the SNN is fed to an ANN, whose role is
to decode ${\bm{y}}$ into the reference signal ${\bm{r}}$. Importantly, we do
not a priori constrain the decoding strategy. In contrast, existing works
typically handcraft target spiking signals ${\bm{y}}$ corresponding to some
fixed decoding rule, such as rate decoding (Fig. 1(a)-(b)). Our main
contribution is a method to train encoding SNN and decoding ANN jointly, with
the goal of discovering automatically efficient time-decoding rules from
spiking signals produced by the encoding SNN.

To optimize SNN encoder and ANN decoder, we consider a data set comprising $N$
pairs $({\bm{x}},{\bm{r}})$ of exogeneous spiking inputs ${\bm{x}}$ and
natural reference signals ${\bm{r}}$ generated from an unknown population
distribution $p({\bm{x}},{\bm{r}})$. The exogeneous inputs are in the form of
arbitrary time vector sequence ${\bm{x}}=({\bm{x}}_{1},\dots,{\bm{x}}_{T})$
with ${\bm{x}}_{t}\in\\{0,1\\}^{N_{X}}$, while the reference signals
${\bm{r}}=({\bm{r}}_{1},\ldots,\bm{r}_{T})$, with
${\bm{r}}_{t}\in\mathbb{R}^{N_{R}}$, define general target signals. In the
image naturalization task illustrated in Fig. 1, the exogenous inputs
${\bm{x}}$ are the spiking signals encoding the input from a DVS camera, with
$N_{X}$ being the number of pixels; and the reference signals ${\bm{r}}$
encode the pixel intensities of the natural image. This can be done for
instant by setting ${\bm{r}}_{t}=\bm{0}_{N_{R}}$ for $t=1,\dots,T-1$, and
${\bm{r}}_{T}$ equal to the grayscale image of $N_{R}$ pixels; or setting
${\bm{r}}_{t}$ to equal the grayscale image of all $t=1,\dots,T$.

We wish to train the encoding SNN to output a spiking representation
${\bm{y}}$ that is maximally informative about the natural target signal
${\bm{r}}$. As we will detail below, the SNN implements a causal stochastic
mapping between latent signal ${\bm{y}}$ and input ${\bm{x}}$ that depends on
a vector of parameters ${\bm{w}}^{e}$ as

$\displaystyle
p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})=\prod_{t=1}^{T}p({\bm{y}}_{t}|{\bm{y}}^{t-1},{\bm{x}}_{t-\tau_{e}}^{t}),$
(3)

where the integer $\tau_{e}$ represents the memory of the encoding SNN and
${\bm{w}}^{e}$ the synaptic weights of the SNN. As a starting point, consider
maximizing the $\tau_{d}$-directed information between ${\bm{y}}$ and
${\bm{r}}$ as

$\displaystyle\max_{{\bm{w}}^{e}}~{}I_{{\bm{w}}^{e}}^{(\tau_{d})}({\bm{y}}\rightarrow{\bm{r}}),$
(4)

for a decoding window $\tau_{d}\geq 0$, over the vector ${\bm{w}}^{e}$ of the
encoding SNN parameters. The use of the $\tau_{d}$-directed information, in
lieu of the mutual information $I({\bm{y}};{\bm{r}})$, reflects the practical
requirement that reference signal ${\bm{r}}_{t}$ be well represented by a
limited window ${\bm{y}}_{t-\tau_{d}}^{t}$ of past outputs of the SNN in a
causal manner. The directed information
$I_{{\bm{w}}^{e}}^{(\tau_{d})}({\bm{y}}\rightarrow{\bm{r}})$ is evaluated with
respect to the marginal $p({\bm{y}},{\bm{r}})$ of the joint distribution

$\displaystyle
p_{{\bm{w}}^{e}}({\bm{x}},{\bm{y}},{\bm{r}})=p({\bm{x}},{\bm{r}})p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}}),$
(5)

combining population distribution and SNN mapping via the chain rule of
probability.

The maximization (4) does not impose any constraint on the complexity of the
spiking representation ${\bm{y}}$. To address this problem, we adopt the
directed information bottleneck (DIB) objective

$R_{\text{DIB}}({\bm{w}}^{e})=I_{{\bm{w}}^{e}}^{(\tau_{d})}({\bm{y}}\rightarrow{\bm{r}})-\beta\cdot
I_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{x}}\rightarrow{\bm{y}}),$ (6)

where $\beta>0$ is a hyperparameter. The penalty term
$I_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{x}}\rightarrow{\bm{y}})$ is the directed
information between the input ${\bm{x}}$ and the spiking representation
${\bm{y}}$, which captures the complexity of the causal encoding done by the
SNN. This term is evaluated with respect to the marginal
$p({\bm{x}},{\bm{y}})=p({\bm{x}})p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})$
of the joint distribution (5), and it reflects the limited memory of the SNN
encoding mapping (3). It can be interpreted as the amount of information about
input ${\bm{x}}$ that is (causally) preserved by the representation
${\bm{y}}$. To the best of our knowledge, the DIB objective is considered in
this paper for the first time.

##  4 Learning to decode through the IB

In order to address the DIB problem (6), we adopt a variational formulation
that relies on the introduction of a decoding network, as in the architecture
illustrated in Fig. 1(d). The decoding network implements the causally
conditional distribution

$\displaystyle
q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})=\prod_{t=1}^{T}q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})$
(7)

between spiking representation ${\bm{y}}$ and natural reference signal
${\bm{r}}$. The decoder (7) is causal, with memory given by integer
$\tau_{d}>1$ and is parametrized by a vector ${\bm{w}}^{d}$. The decoder can
be implemented in several ways. A simple approach, that we follow in the
remainder of this paper, is to use a model
$q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})=q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{y}}_{t-\tau_{d}}^{t})$,
where $q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{y}}_{t-\tau_{d}}^{t})$ is an ANN
with inputs given by a window ${\bm{y}}_{t-\tau_{d}}^{t}$ of $\tau_{d}$
samples from the spiking representation ${\bm{y}}$. Alternatively, one could
use an RNN to directly model the kernel
$q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})$.
Note that the introduction of the decoder network (7) is consistent with the
use of the $\tau_{d}$-directed information in (6). We emphasize that the use
of an ANN, or RNN, for decoding, is essential in order to allow the reference
sequence ${\bm{r}}$ to be a natural signal.

With such a network, using a standard variational inequality [9], we bound the
two terms in (6) to obtain the variational DIB (VDIB) loss
$\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})$ as an upper bound on
the negative DIB objective (6)

$\displaystyle\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})=$
$\displaystyle\underbrace{\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}},{\bm{r}})}\Bigg{[}-\log
q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})\Bigg{]}}_{\text{average
decoder's log-
loss}}+\beta\cdot\mathrm{E}_{p({\bm{x}})}\Bigg{[}\underbrace{\text{KL}(p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})||q({\bm{y}}))}_{\text{information
theoretic regularization}}\Bigg{]},$ (8)

in which we have introduced an arbitrary “prior” distribution
$q({\bm{y}})=\prod_{t}q({\bm{y}}_{t}|{\bm{y}}^{t-1})$ on the spiking
representation. As in the implementation in [9] for the standard IB criterion,
we will consider $q({\bm{y}})$ to be fixed, although it can potentially be
optimized on. To enforce sparsity of the encoder’s outputs, the “prior”
$q({\bm{y}})$ can be chosen as a sparse Bernoulli distribution [4], i.e.,
$q({\bm{y}})=\prod_{t=1}^{T}\text{Bern}({\bm{y}}_{t}|p)$ for some small
probability $p$. The derivation of the bound can be found in Appendix A.1.

The learning criterion (8) enables the joint training of encoding SNN and
decoding network via the problem

$\displaystyle\min_{{\bm{w}}^{e},{\bm{w}}^{d}}\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d}),$
(9)

which we address via SGD. To this end, we now give general expressions for the
gradients of the VDIB criterion (8) with respect to the encoder and decoder
weights by leveraging the REINFORCE Monte Carlo gradients. In the next
section, we will then detail how the gradient with respect to the encoding SNN
weights ${\bm{w}}^{e}$ can be approximated by integrating SG [3] with the
probabilistic approach [4]. To start, the VDIB loss in (8) can equivalently be
stated as

$\displaystyle\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})=\mathrm{E}_{\underbrace{p({\bm{x}},{\bm{r}})}_{\text{population
distribution}}}\mathrm{E}_{\underbrace{p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})}_{\text{enc.
SNN}}}\Bigg{[}\underbrace{-\log
q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})}_{:=~{}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})}+\beta\cdot\underbrace{\log\Bigg{(}\frac{p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})}{q({\bm{y}})}}_{:=~{}\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})}\Bigg{)}\Bigg{]},$
(10)

where we have defined the decoder loss
$\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})$ and the encoder loss
$\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})$, as shown in (10). The gradient with
respect to the encoder weights ${\bm{w}}^{e}$ can be obtained via the
REINFORCE gradient technique [22, Ch. 8][23] as

$\displaystyle\nabla_{{\bm{w}}^{e}}\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})=\mathrm{E}_{p({\bm{x}},{\bm{r}})}\mathrm{E}_{p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})}\bigg{[}\Big{(}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})+\beta\cdot\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})\Big{)}\cdot\nabla_{{\bm{w}}^{e}}\log
p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})\bigg{]},$ (11)

while the gradient with respect to the decoder weights ${\bm{w}}^{d}$ can be
directly computed as

$\displaystyle\nabla_{{\bm{w}}^{d}}\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})=\mathrm{E}_{p({\bm{x}},{\bm{r}})}\mathrm{E}_{p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})}\bigg{[}\nabla_{{\bm{w}}^{d}}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})\bigg{]}.$
(12)

Unbiased estimates of the gradients (11)-(12) can be evaluated by using one
sample $({\bm{x}},{\bm{r}})\sim p({\bm{x}},{\bm{r}})$ from the data set, along
with a randomly generated SNN output ${\bm{y}}\sim
p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})$. We obtain

$\displaystyle\nabla_{{\bm{w}}^{e}}\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})\approx\Big{(}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})+\beta\cdot\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})\Big{)}\cdot\nabla_{{\bm{w}}^{e}}\log
p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})$ (13)
$\displaystyle\nabla_{{\bm{w}}^{d}}\mathcal{L}_{\text{VDIB}}({\bm{w}}^{e},{\bm{w}}^{d})\approx\nabla_{{\bm{w}}^{d}}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}}).$
(14)

While the gradient
$\nabla_{{\bm{w}}^{d}}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})$ can be computed
using standard backpropagation on the decoding ANN, the gradient
$\nabla_{{\bm{w}}^{e}}\log p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})$
for the encoder weights depends on the SNN model, and is detailed in the next
section.

As we will describe in the next section, and as illustrated in Fig. 1(d), the
learning rule for the SNN involves the common learning signal
$\big{(}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})+\beta\cdot\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})\big{)}$
which includes two feedback signals, one from the ANN decoder – the loss
$\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})$ – and one from the SNN encoder,
namely $\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})$. The latter imposes a penalty
for the distribution of the representation ${\bm{y}}$ that depends on its
“distance”, in terms of log-likelihood ratio, with respect to the reference
distribution $q({\bm{y}})$. Apart from the described common learning signals,
which is a hallmark of probabilistic learning rules [4, 24], the update also
involves per-neuron terms derived as in SG-based rules. The proposed algorithm
is detailed in Appendix A.1 and derived next.

##  5 SNN Model

The encoding SNN is a directed network $\mathcal{N}$ of spiking neurons for
which at each time-step $t=1,\dots,T$, each neuron $i\in\mathcal{N}$ outputs a
binary value $z_{i,t}\in\\{0,1\\}$. The set of neurons can be partitioned as
$\mathcal{N}=\mathcal{H}\cup\mathcal{Y}$, where $\mathcal{H}$ denotes the set
of $N_{H}$ hidden neurons, and $\mathcal{Y}$ the set of $N_{Y}$ readout
neurons that output spiking signal ${\bm{y}}$ over $T$ time instants. The SNN
defines the stochastic mapping
$p^{(\tau_{e})}_{{\bm{w}}^{e}}({\bm{y}}||{\bm{x}})$ in (3) between exogenous
input sequence ${\bm{x}}$ and the output ${\bm{y}}$ of the readout neurons.
For each neuron $i\in\mathcal{N}$, we denote as $\mathcal{P}_{i}$ the set of
pre-synaptic neurons, i.e., the set of neurons having directed synaptic links
to neuron $i$. With a slight abuse of notation, we include the $N_{X}$
exogenous inputs in this set. In order to enable the use of information-
theoretic criteria via the stochastic mapping (3), as well as the efficiency
of SG-based training, we model the readout neurons in $\mathcal{Y}$ using the
(probabilistic) Spike Response Model (SRM) [25] with stochastic threshold [4],
while the hidden neurons in $\mathcal{H}$ follow the standard (deterministic)
SRM. Accordingly, the spiking mechanism of every neuron $i\in\mathcal{N}$ in
the SNN depends on the spiking history of pre-synaptic neurons and of the
post-synaptic neuron itself through the neuron’s membrane potential [25, 4,
24]

$\displaystyle
u_{i,t}=\sum_{j\in\mathcal{P}_{i}}w^{e}_{ij}\big{(}\alpha_{t}\ast
z_{j,t}\big{)}+w^{e}_{i}\big{(}\beta_{t}\ast
z_{i,t-1}\big{)}+\bar{w}_{i}^{e}.$ (15)

In (15), the contribution of each pre-synaptic neuron $j$ is given by the
synaptic trace $\alpha_{t}\ast z_{j,t}$, where the sequence $\alpha_{t}$ is
the synaptic spike response and $z_{j,t}$ is the pre-synaptic neuron’s output,
multiplied by the synaptic weight $w^{e}_{ij}$. The membrane potential is also
affected by the post-synaptic trace $\beta_{t}\ast z_{i,t-1}$, where
$\beta_{t}$ is the feedback response and $z_{i,t}$ is the post-synaptic
neuron’s output, which is weighted by $w^{e}_{i}$. The spike response
$\alpha_{t}$ and feedback response $\beta_{t}$ can be different for visible
and hidden neurons (see Appendix A.2), and they have memory limited for
$\tau_{e}$ samples, i.e., $\alpha_{t}=\beta_{t}=0$ for $t>\tau_{e}$. Finally,
$\bar{w}_{i}^{e}$ is a learnable bias parameter. Choices for these functions
and details on the neuron models are provided in Appendix A.2.

#### Deterministic SRM.

For each neuron $i\in\mathcal{N}$, the binary output $z_{i,t}\in\\{0,1\\}$ at
time $t$ depends on its membrane potential $u_{i,t}$ as

$\displaystyle z_{i,t}=\Theta(u_{i,t}),$ (16)

where $\Theta(\cdot)$ is the Heaviside step function: A spike is emitted when
the membrane potential $u_{i,t}$ is positive.

#### SRM with stochastic threshold.

For each neuron $i\in\mathcal{Y}$ in the readout layer, the spiking
probability, when conditioned on the overall spiking history, depends on the
membrane potential as

$\displaystyle p_{{\bm{w}}^{e}_{i}}(y_{i,t}=1|u_{i,t})=\sigma(u_{i,t}),$ (17)

where $\sigma(x)=1/(1+\exp{(-x)})$ is the sigmoid function, and
${\bm{w}}^{e}_{i}=\\{w^{e}_{ij},w^{e}_{i},\bar{w}_{i}^{e}\\}$ is the per-
neuron set of parameters. This yields the encoding stochastic mapping (3) as

$\displaystyle
p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})=\prod_{i\in\mathcal{Y}}\prod_{t=1}^{T}p_{{\bm{w}}^{e}_{i}}(y_{i,t}||u_{i,t}).$
(18)

#### Gradients.

From Sec. 4, given (13)-(14), we need to compute the gradients
$\nabla_{{\bm{w}}^{e}}\log p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})$.
For the readout neurons, this can be directly done as [4]

$\displaystyle\nabla_{w^{e}_{ij}}\log
p_{{\bm{w}}_{i}^{e}}(y_{i,t}=1|u_{i,t})=\big{(}\alpha_{t}\ast
z_{j,t}\big{)}\big{(}y_{i,t}-\sigma(u_{i,t})\big{)}=e_{ij,t},$ (19)

where we recall that $z_{j,t}$ is the spiking signal of pre-synaptic neuron
$j$. In (19), we have defined the pre-synaptic eligibility trace $e_{ij,t}$.
Learning rules can similarly derived for weights $w^{e}_{i}$ and
$\bar{w}^{e}_{i}$, with corresponding eligibility traces $e_{i,t}$ and
$\bar{e}_{i,t}$, and we define per-neuron eligibility trace
${\bm{e}}_{i,t}=\\{e_{ij,t},e_{i,t},\bar{e}_{i,t}\\}$. For neurons
$i\in\mathcal{H}$ in the hidden layer, we rely on e-prop [3] – an SG-based
technique – to obtain the approximation

$\displaystyle\nabla_{w^{e}_{ij}}\log
p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})\approx\sum_{t=1}^{T}L_{i,t}e_{ij,t},$
(20)

with the learning signal

$\displaystyle
L_{i,t}=\sum_{k\in\mathcal{Y}}B_{ik}\big{(}y_{k,t}-\sigma(u_{k,t})\big{)}$
(21)

characterized by random, but fixed, weights $B_{ik}$ following direct feedback
alignment [26], and eligibility trace

$\displaystyle
e_{ij,t}=\Theta^{\prime}(u_{i,t}-\vartheta)\big{(}\alpha_{t}\ast
z_{j,t}\big{)}.$ (22)

The approximation (20) comes from: (i) truncating the BPTT by ignoring the
dependence of the log-loss $-\log p_{{\bm{w}}_{i}^{e}}(y_{i,t}|u_{i,t})$ at
time $t$ on the weight $w_{i}^{e}$ due to the previous time instants
$t^{\prime}<t$; (ii) replacing backpropagation through neurons, via feedback
alignment; and (iii) using the derivative of a surrogate function
$\Theta^{\prime}(\cdot)$ [12]. For example, with a sigmoid surrogate gradient,
we have $\Theta^{\prime}(\cdot)=\sigma^{\prime}(\cdot)$.

The overall algorithm, referred to as VDIB, is summarized in Appendix A.1. As
anticipated in the previous sections, the updates (31) and (32) contain
elements from both probabilistic and SG-based learning rules. In particular,
the global learning signal
$\big{(}\ell_{{\bm{w}}^{d}}({\bm{y}},{\bm{r}})+\beta\cdot\ell_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})\big{)}$
relates to the global signal used in probabilistic learning [4, 24], whilst
per-neuron error signal (21) is a feature of SG schemes [12, 3].

##  6 Experiments

We now evaluate the proposed VDIB solution in an experimental setting. We
specifically focus on two tasks: predictive coding and naturalization of
spiking signals. The former task highlights the aspect of causal, online,
encoding and decoding, and relates closely with work in theoretical
neuroscience [21]. In contrast, the latter task emphasizes the role of the
proposed approach for neuromorphic computing [6]. Specifically, the
naturalization of spiking inputs, e.g., obtained from a DVS camera [5],
enables the integration of neuromorphic sensors with conventional digital
sensors and processors [6].

As the proposed hybrid architecture gives flexibility in the choice of the
decoding ANN, we have considered three options: (i) a simple logistic
regression model, (ii) a multilayer perceptron (MLP) with one hidden layer,
and (iii) a convolutional neural network (CNN) with two convolutional layers
and one fully connected layer. As discussed in Sec. 3, at each time $t$, the
decoding ANN is given the sequence ${\bm{y}}_{t-\tau_{d}}^{t}$ of $\tau_{d}$
time samples produced by the encoding SNN. As a benchmark, we also consider
conventional rate decoding, whereby time samples within the sequence
${\bm{y}}_{t-\tau_{d}}^{t}$ are summed over time, yielding
$\sum_{t^{\prime}=t-\tau_{d}}^{t}{\bm{y}}_{t^{\prime}}$, before being fed to
the ANN.

By considering encoding networks with a layered architecture, our
implementation of both ANN and SNN makes use of automatic differentiation and
general-purpose deep learning libraries [27] and Tesla V100 GPUs.

###  6.1 Efficient predictive coding

Figure 2: Left: Illustration of the exogenous inputs (top), and accuracy as a
function of the target lag $\delta$ (bottom). Right: Encoding network spike
rates with lag $\delta=-2$ versus regularization constant $\beta$ (top), and
corresponding accuracy (bottom). Shaded areas represent standard deviations.

In this first experiment, we explore the use of the proposed hybrid
variational autoencoder for a predictive coding problem. Following the
experiment proposed in [21], we consider exogenous signals consisting of $20$
correlated spiking signals that represent two independent auto-regressive
“drifting Gaussian blobs”. At each time $t$, the position
$p_{t}\in\\{1,\dots,20\\}$ of each “blob” across the exogenous inputs follows
a (wrapped) Gaussian path $\mathcal{N}(\theta_{t},\sigma)$, with $\theta_{t}$
evolving according to an AR(2) process $\theta_{t}=\theta_{t-1}+v_{t}$, with
$v_{t}=av_{t-1}-b\epsilon_{t}$, $\epsilon_{t}\sim\mathcal{N}(0,1)$,
$\sigma=0.45$, $a=0.9$ and $b=0.14$. A realization of the exogenous inputs can
be found in Fig. 2(left, top). The reference signal ${\bm{r}}_{t}$ at time-
step $t$ consists of the one-hot encoding of the positions of the two blobs at
time $t+\delta$ for some time lag $\delta\in\mathbb{Z}$. A one layer encoder
compresses the $20$ exogenous spiking signals into $10$ spiking signals, which
are decoded by a softmax regressor (playing the role of the decoding ANN). We
choose the reference distribution $q({\bm{y}})$ to be i.i.d.
$\text{Bern}(0.2)$ to enforce temporal sparsity. We set $\tau_{e}=5$, and
$\tau_{d}=5$ and, where not mentioned otherwise, $\beta=1$. Training is done
over $50,000$ randomly generated examples, of length $T=100$. Testing is done
on a single example of length $T=1,000$. Results are presented in terms of
mean squared error (MSE), and averaged over $5$ trials with random
initializations. In the right panels, it can be seen that by modulating the
value of $\beta$, one can explore the trade-off between temporal sparsity and
accuracy. Overall, the approach is seen to yield compact spiking
representations that enable the reconstruction of the input signal via online
decoding. Finally, we provide a study of the effect of the window sizes
$\tau_{e}$ and $\tau_{d}$ in Appendix A.3.

###  6.2 Naturalizing MNIST images

We now consider image naturalization. To start, we obtain exogenous input
spiking signals ${\bm{x}}$ by encoding MNIST images using either Poisson
encoding (see [25]) or time-to-first-spike encoding (as in [7]). The relevance
signal is defined as ${\bm{r}}_{t}=\bm{0}_{N_{X}}$ for $t=1,\dots,T-1$, and
${\bm{r}}_{T}$ is the original image. Results are provided in terms of the
accuracy of the convolutional LeNet ANN trained on the original MNIST dataset
[28]. Specifically, we feed the decoded natural images to the LeNet ANN, and
report its accuracy. The hybrid variational autoencoder compresses the $784$
exogenous spiking signals into $256$ spiking signals at the readout layer of
the SNN. We set $\tau_{e}=\tau_{d}=T=30$, and $\beta=0.001$. Training is
carried over $200,000$ examples, and testing is done on the $10,000$ images
from the test dataset. Other experimental details can be found in Appendix
A.3.

Figure 3: Top: Exogenous spiking inputs from the MNIST-DVS dataset summed over time samples (first row), and reconstructed natural images via VDIB (second row). Bottom: T-SNE representation of the spiking signals produced by the SNN with rate decoding (left), and with VDIB (right). Table 1: Comparison of classification accuracies for an MLP ANN decoder under various encoding strategies (for the MNIST digits at the input of the SNN) and decoding schemes (for the output of the SNN). Encoding \ Decoding | Rate | VDIB
---|---|---
Rate | $60.50\pm 0.32\%$ | $\bm{91.32\pm 0.27\%}$
Time-to-First-Spike | $86.43\pm 0.84\%$ | $\bm{92.82\pm 0.25\%}$

We explore the impact of various encoding and decoding strategies in Table 1,
when using a decoding MLP. From Table 1, we note that, although rate decoding
provides clear benefits in reducing the size of inputs to the decoder, this is
at the cost of a large drop in accuracy. The clear performance advantages for
time decoding validate our working hypothesis regarding the importance to move
away from handcrafted decoding strategies.

###  6.3 Naturalizing MNIST-DVS images

Finally, we evaluate VDIB on the MNIST-DVS dataset [29]. Following [30, 31,
32, 16], we crop images to $26\times 26$ pixels, and to $2$ seconds. We also
use a sampling period of $10$ ms, which yields $T=200$, and evaluate the
impact of coarser sampling rates on the proposed method. To obtain natural
images, the reference signal ${\bm{r}}$ used during training is given by a
single image of the corresponding class from the MNIST dataset. In this
experiment, the $676$ exogenous signals from MNIST-DVS images are compressed
by the encoder into $256$ spiking signals, which are fed to a convolutional
decoder. We set $\tau_{e}=\tau_{d}=T$, and $\beta=0.001$. Training is carried
over $100,000$ examples, and testing is done on the $10,000$ images from the
test dataset.

Fig. 3(top row) illustrates with some examples that VDIB can effectively
naturalize MNIST-DVS images. In the bottom row, we show T-SNE [33]
representations of encoded spiking signals produced by the SNN when using
VDIB, and compare them with those produced when considering rate decoding. The
latter is seen to be unable to produce representations from which class
patterns can be identified, since all the representations are superimposed.

##  7 Conclusion

In this work, we have introduced VDIB, a variational learning algorithm for a
hybrid SNN-ANN autoencoder that is based on a directed, causal, variant of the
information bottleneck. The proposed system enables a shift from conventional
handcrafted decoding rules to optimized learned decoding. VDIB can find
applications in neuromorphic sensing through naturalization of spiking
signals, and generally as a broad framework for the design of training
algorithms for SNNs. This includes applications in mobile edge computing,
offering opportunities in, e.g., personal healthcare. Limitations of this work
include the choice of ANNs for decoders, where the usage of RNNs may result in
improvements. Another line for future work consists in learning the reference
distribution $q({\bm{y}})$ instead of fixing it.

## Acknowledgments and Disclosure of Funding

This work was supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 725731) and by Intel Labs via the Intel Neuromorphic Research
Community.

The authors acknowledge use of the Joint Academic Data science Endeavour
(JADE) HPC Facility (http://www.jade.ac.uk/), and use of the research
computing facility at King’s College London, Rosalind
(https://rosalind.kcl.ac.uk).

## References

  * [1] Mark Humphries,  The Spike: An Epic Journey Through the Brain in 2.1 Seconds,  Princeton University Press, 2021.
  * [2] Mike Davies, Andreas Wild, Garrick Orchard, Yulia Sandamirskaya, Gabriel A. Fonseca Guerra, Prasad Joshi, Philipp Plank, and Sumedh R. Risbud,  “Advancing Neuromorphic Computing with Loihi: A Survey of Results and Outlook,”  Proceedings of the IEEE, vol. 109, no. 5, pp. 911–934, 2021.
  * [3] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Robert Salaj, Darjan Legenstein, and Wolfgang Maass,  “A Solution to the Learning Dilemma for Recurrent Networks of Spiking Neurons,”  Nature Communications, vol. 11, no. 3625, 2020.
  * [4] Hyeryung Jang, Osvaldo Simeone, Brian Gardner, and André Grüning,  “An Introduction to Probabilistic Spiking Neural Networks: Probabilistic Models, Learning Rules, and Applications,”  IEEE Sig. Proc. Magazine, vol. 36, no. 6, pp. 64–77, 2019.
  * [5] Patrick Lichtsteiner, Christoph Posch, and Tobi Delbruck,  “A 128 x 128 120db 30mw Asynchronous Vision Sensor that Responds to Relative Intensity Change,”  in Proc. Int. Solid State Circuits Conf., 2006, pp. 2060–2069.
  * [6] Dennis Robey, Wesley Thio, Herbert H. C. Iu, and Jason Kamran Eshraghian,  “Naturalizing Neuromorphic Vision Event Streams Using GANs,”  CoRR, vol. abs/2102.07243, 2021.
  * [7] Alireza Bagheri, Osvaldo Simeone, and Bipin Rajendran,  “Training Probabilistic Spiking Neural Networks with First-To-Spike Decoding,”  in 2018 IEEE Int. Conf. on Acoustics, Speech and Sig. Processing (ICASSP), 2018, pp. 2986–2990.
  * [8] Naftali Tishby, Fernando C. Pereira, and William Bialek,  “The Information Bottleneck Method,”  in The 37th annual Allerton Conf. on Communication, Control, and Computing, 1999, pp. 368–377.
  * [9] Alex Alemi, Ian Fischer, Josh Dillon, and Kevin Murphy,  “Deep Variational Information Bottleneck,”  in ICLR, 2017.
  * [10] Gerhard Kramer,  Directed Information for Channels with Feedback,  Hartung-Gorre, 1998.
  * [11] Emre Neftci, Charles Augustine, Somnath Paul, and Georgios Detorakis,  “Event-Driven Random Back-propagation: Enabling Neuromorphic Deep Learning Machines,”  Frontiers in neuroscience, vol. 11, pp. 324, 2017.
  * [12] Emre Neftci, Hesham Mostafa, and Friedemann Zenke,  “Surrogate Gradient Learning in Spiking Neural Networks: Bringing the Power of Gradient-Based Optimization to Spiking Neural Networks,”  IEEE Sig. Proc. Magazine, vol. 36, no. 6, pp. 51–63, 2019.
  * [13] Jacques Kaiser, Hesham Mostafa, and Emre Neftci,  “Synaptic plasticity dynamics for deep continuous local learning (DECOLLE),”  Frontiers in Neuroscience, vol. 14, pp. 424, 2020.
  * [14] Sumit Bam Shrestha and Garrick Orchard,  “SLAYER: Spike Layer Error Reassignment in Time,”  in Advances in Neural Information Processing Systems, 2018, vol. 31.
  * [15] Friedemann Zenke and Surya Ganguli,  “SuperSpike: Supervised Learning in Multilayer Spiking Neural Networks,”  Neural Computation, vol. 30, no. 6, pp. 1514–1541, 2018.
  * [16] Hyeryung Jang, Nicolas Skatchkovsky, and Osvaldo Simeone,  “VOWEL: A Local Online Learning Rule for Recurrent Networks of Probabilistic Spiking Winner-Take-All Circuits,”  arXiv preprint arXiv:2004.09416, 2020.
  * [17] Kenneth Stewart, Andreea Danielescu, Lazar Supic, Timothy Shea, and Emre Neftci,  “Gesture Similarity Analysis on Event Data Using a Hybrid Guided Variational Auto Encoder,”  arXiv:2104.00165, 2021.
  * [18] Naftali Tishby and Noga Zaslavsky,  “Deep Learning and the Information Bottleneck Principle,”  in 2015 IEEE Information Theory Workshop (ITW), 2015, pp. 1–5.
  * [19] Stefan Klampfl, Robert Albin Legenstein, and Wolfgang Maass,  “Spiking Neurons Can Learn to Solve Information Bottleneck Problems and to Extract Independent Components,”  Neural computation, vol. 21, no. 4, pp. 911–959, 2009.
  * [20] Lars Buesing and Wolfgang Maass,  “A Spiking Neuron as Information Bottleneck,”  Neural Computation, vol. 22, no. 8, pp. 1961–1992, 08 2010.
  * [21] Matthew Chalk, Olivier Marre, and Gašper Tkačik,  “Toward a Unified Theory of Efficient, Predictive, and Sparse Coding,”  Proceedings of the National Academy of Sciences, vol. 115, no. 1, pp. 186–191, 2018.
  * [22] Osvaldo Simeone,  “A Brief Introduction to Machine Learning for Engineers,”  Foundations and Trends® in Signal Processing, vol. 12, no. 3-4, pp. 200–431, 2018.
  * [23] Andriy Mnih and Danilo Rezende,  “Variational Inference for Monte Carlo Objectives,”  in International Conference on Machine Learning, 2016, pp. 2188–2196.
  * [24] Hyeryung Jang, Nicolas Skatchkovsky, and Osvaldo Simeone,  “Spiking Neural Networks – Parts I, II, and III,”  https://arxiv.org/abs/2010.14208, 2020.
  * [25] Wulfram Gerstner and Werner M Kistler,  Spiking Neuron Models: Single Neurons, Populations, Plasticity,  Cambridge University Press, 2002.
  * [26] Arild Nøkland and Lars Hiller Eidnes,  “Training Neural Networks with Local Error Signals,”  in Proceedings of the 36th International Conference on Machine Learning. 2019, vol. 97, pp. 4839–4850, PMLR.
  * [27] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala,  “Pytorch: An Imperative Style, High-Performance Deep Learning Library,”  in Advances in Neural Information Processing Systems 32, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., pp. 8024–8035. Curran Associates, Inc., 2019.
  * [28] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel,  “Backpropagation Applied to Handwritten Zip Code Recognition,”  Neural Computation, vol. 1, no. 4, pp. 541–551, 1989.
  * [29] Teresa Serrano-Gotarredona and Bernabé Linares-Barranco,  “Poker-DVS and MNIST-DVS. Their History, How They Were Made, and Other Details,”  Frontiers in neuroscience, vol. 9, pp. 481, 2015.
  * [30] Bo Zhao, Ruoxi Ding, Shoushun Chen, Bernabe Linares-Barranco, and Huajin Tang,  “Feedforward Categorization on AER Motion Events Using Cortex-Like Features in a Spiking Neural Network,”  IEEE Trans. on Neural Net. and Learning Systems, vol. 26, no. 9, pp. 1963–1978, 2014.
  * [31] James A Henderson, TingTing A Gibson, and Janet Wiles,  “Spike event based learning in neural networks,”  arXiv preprint arXiv:1502.05777, 2015.
  * [32] N. Skatchkovsky, H. Jang, and O. Simeone,  “Federated neuromorphic learning of spiking neural networks for low-power edge intelligence,”  in Proc. of International Conference on Acoustics, Speech and Signal Processing. IEEE, 2020.
  * [33] Laurens van der Maaten and Geoffrey Hinton,  “Visualizing Data using t-SNE,”  Journal of Machine Learning Research, vol. 9, no. 86, pp. 2579–2605, 2008.
  * [34] Jonathan W. Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M. Litke, E.J. Chichilnisky, and Eero P. Simoncelli,  “Spatio-Temporal Correlations and Visual Signalling in a Complete Neuronal Population,”  Nature, vol. 454, no. 7207, pp. 995–999, 2008.

##  Appendix A Appendix

###  A.1 Derivation of the bound

As described in Sec 4, using a decoding network defined by the causally
conditional distribution $q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})$
from eq. (7), we use a standard variational inequality [9] to lower bound the
mutual information at each step $t=1,\dots,T$

$\displaystyle
I_{{\bm{w}}^{e}}^{(\tau_{d})}({\bm{y}}_{t-\tau_{d}}^{t};{\bm{r}}_{t}|{\bm{r}}^{t-1})$
$\displaystyle\geq\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}}^{t},{\bm{r}}^{t})}\Bigg{[}\log\frac{q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})}{p({\bm{r}}_{t}|{\bm{r}}^{t-1})}\Bigg{]}$
(23)
$\displaystyle=\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}}^{t},{\bm{r}}^{t})}\Bigg{[}\log
q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})\Bigg{]}+H({\bm{r}}_{t}|{\bm{r}}^{t-1}).$
(24)

Summing over the time-steps, we then obtain

$\displaystyle I_{{\bm{w}}^{e}}^{(\tau_{d})}({\bm{y}}\rightarrow{\bm{r}})$
$\displaystyle\geq\sum_{t=1}^{T}\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}}^{t},{\bm{r}}^{t})}\Bigg{[}\log
q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})\Bigg{]}+\sum_{t=1}^{T}H({\bm{r}}_{t}|{\bm{r}}^{t-1})$
(25)
$\displaystyle=\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}},{\bm{r}})}\Bigg{[}\sum_{t=1}^{T}\log
q_{{\bm{w}}^{d}}({\bm{r}}_{t}|{\bm{r}}^{t-1},{\bm{y}}_{t-\tau_{d}}^{t})\Bigg{]}+H({\bm{r}})$
(26)
$\displaystyle=\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}},{\bm{r}})}\Bigg{[}\log
q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})\Bigg{]}+H({\bm{r}}),$ (27)

where the entropy of the target signals $H({\bm{r}})$ is independent of the
model parameters ${\bm{w}}^{e}$ and can be ignored for the purpose of
training.

To bound the second information term in (6), we introduce an arbitrary
auxiliary distribution $q({\bm{y}})=\prod_{t}q({\bm{y}}_{t}|{\bm{y}}^{t-1})$
on the readout layer. As in the implementation in [9], we will consider it to
be fixed, although it can potentially be optimized on. Using a standard
variational inequality [9], we then obtain the upper bound the mutual
information of the encoder for each $t=1,\dots,T$ as

$\displaystyle I_{{\bm{w}}^{e}}({\bm{x}}^{t};{\bm{y}}^{t})$
$\displaystyle\leq\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}}^{t},{\bm{x}}^{t})}\Bigg{[}\log\frac{p_{{\bm{w}}^{e}}({\bm{y}}_{t}|{\bm{x}}^{t}_{t-\tau_{e}},{\bm{y}}^{t-1})}{q({\bm{y}}_{t}|{\bm{y}}^{t-1})}\Bigg{]}.$
(28)

Summing over time-steps $t=1,\dots,T$, we obtain

$\displaystyle I_{{\bm{w}}^{e}}({\bm{x}}\rightarrow{\bm{y}})$
$\displaystyle\leq\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}},{\bm{x}})}\Bigg{[}\log\frac{p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})}{q({\bm{y}})}\Bigg{]}.$
(29)

Combining bounds (25)-(29), we obtain the variational DIB loss
$\mathcal{L}_{VDIB}({\bm{w}}^{e},{\bm{w}}^{d})$ as an upper bound on the
negative DIB objective

$\displaystyle\mathcal{L}_{VDIB}({\bm{w}}^{e},{\bm{w}}^{d})$
$\displaystyle=\mathrm{E}_{p_{{\bm{w}}^{e}}({\bm{y}},{\bm{r}})}\Bigg{[}-\log
q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}||{\bm{y}})\Bigg{]}+\beta\cdot\mathrm{E}_{p({\bm{x}})}\Bigg{[}\text{KL}(p_{{\bm{w}}^{e}}^{(\tau_{e})}({\bm{y}}||{\bm{x}})||q({\bm{y}}))\Bigg{]}.$
(30)

This completes the derivation. The proposed method is detailed in Algorithm 1.

Input: Exogeneous spiking signal $\mathbf{x}$, reference natural signal
$\mathbf{r}$, reference distribution $q({\bm{y}})$, learning rate $\eta$,
moving average parameter $\kappa$, regularization strength $\beta$

Output: Learned model parameters ${\bm{w}}^{e},{\bm{w}}^{d}$

1 initialize parameters ${\bm{w}}^{e},{\bm{w}}^{d}$;

2 for  _each iteration_ do

3 draw a sample $({\bm{x}},{\bm{r}})$ from the data set

4 for  _each time $t=1,\dots,T$_ do

5 \- generate spike outputs from the encoding SNN ${\bm{y}}_{t}\sim
p_{{\bm{w}}^{e}}({\bm{y}}_{t}||{\bm{x}}^{t})$

6 \- a central processor collects the log-probabilities $\log
p_{{\bm{w}}^{e}}(y_{i,t}||{\bm{x}}^{t})$ for all readout neurons
$i\in\mathcal{Y}$ and updates the encoder loss
$\displaystyle\ell_{{\bm{w}}^{e}}({\bm{y}}_{t},{\bm{x}}_{t},{\bm{h}}_{t})=\sum_{i\in\mathcal{Y}}\log\Bigg{(}\frac{p_{{\bm{w}}^{e}}(y_{i,t}||{\bm{x}}^{t})}{q(y_{i,t})}\Bigg{)}$

7 \- generate outputs from the decoding network ${\bm{r}}_{t}\sim
q_{{\bm{w}}^{d}}({\bm{r}}_{t}||{\bm{y}}^{t})$

8 \- a central processor collects the log-probabilities $\log
q_{{\bm{w}}^{d}}({\bm{r}}||{\bm{y}}^{t})$ and updates the decoder loss
$\displaystyle\ell_{{\bm{w}}^{d}}({\bm{y}}_{t},{\bm{r}}_{t})=-\log
q_{{\bm{w}}^{d}}({\bm{r}}_{t}||{\bm{y}}^{t})$

9 for  _each encoding network readout neuron $i\in\mathcal{Y}$_ do

10 compute update $\Delta_{i,t}$ using (19) as
$\displaystyle\Delta_{i,t}={\bm{e}}_{i,t}$

11 end for

12 for  _each encoding network hidden neuron $i\in\mathcal{H}$_ do

13 compute update $\Delta_{i,t}$ using (20) as
$\displaystyle\Delta_{i,t}=L_{i,t}\cdot{\bm{e}}_{i,t}$

14 end for

15 update the encoding SNN model parameters as
$\displaystyle{\bm{w}}^{e}\leftarrow{\bm{w}}^{e}-\eta\cdot\big{(}\ell_{{\bm{w}}^{d}}({\bm{y}}_{t},{\bm{r}}_{t})+\beta\cdot\ell_{{\bm{w}}^{e}}({\bm{y}}_{t},{\bm{x}}_{t},{\bm{h}}_{t})\big{)}\cdot\Delta_{t}$
(31)

16 update decoding ANN model parameters as
$\displaystyle{\bm{w}}^{d}\leftarrow{\bm{w}}^{d}-\eta\cdot\nabla_{{\bm{w}}^{d}}\ell_{{\bm{w}}^{d}}({\bm{y}}_{t},{\bm{r}}_{t})$
(32)

17 end for

18

19 end for

Algorithm 1 Learn to Time-Decode via Variational Directed IB (VDIB)

###  A.2 Neuron models

#### Deterministic SRM.

For hidden neurons $i\in\mathcal{H}$, synaptic and feedback filters are chosen
as the alpha-function spike response
$\alpha_{t}=\exp(-t/\tau_{\text{mem}})-\exp(-t/\tau_{\text{syn}})$ and the
exponentially decaying feedback filter $\beta_{t}=-\exp(-t/\tau_{\text{ref}})$
for $t\geq 1$ with some positive constants
$\tau_{\text{mem}},\tau_{\text{syn}}$, and $\tau_{\text{ref}}$. It is then
common to compute the membrane potential (15) using the following set of
autoregressive equations [24]

$\displaystyle p_{j,t}$
$\displaystyle=\exp\big{(}-1/\tau_{\text{mem}}\big{)}p_{j,t-1}+q_{j,t-1},$
(33a) $\displaystyle\text{with}~{}~{}q_{j,t}$
$\displaystyle=\exp\big{(}-1/\tau_{\text{syn}}\big{)}q_{j,t-1}+s_{j,t-1},$
(33b) $\displaystyle\text{and}~{}~{}r_{i,t}$
$\displaystyle=\exp\big{(}-1/\tau_{\text{ref}}\big{)}r_{i,t-1}+s_{i,t-1}.$
(33c)

Equations (33a)-(33b) define the computation of the pre-synaptic contribution
via a second-order autoregressive (AR) filter, while the feedback of neuron
$i$ is computed in (33c) via a first-order AR filter.

#### SRM with stochastic threshold.

For readout neurons, we associate each synapse with $K_{a}\geq 1$ spike
responses $\\{\alpha_{t}^{k}\\}_{k=1}^{K_{a}}$ and corresponding weights
$\\{w^{e}_{ij,k}\\}_{k=1}^{K_{a}}$, so that the contribution from pre-synaptic
neuron $j$ to the membrane potential in (15) can be written as the double sum
[34]

$\displaystyle\sum_{j\in\mathcal{P}_{i}}\sum_{k=1}^{K_{a}}w^{e}_{ij,k}\big{(}\alpha_{t}^{k}\ast
z_{j,t}\big{)}.$ (34)

By choosing the spike responses $\alpha_{t}^{k}$ to be sensitive to different
synaptic delays, this parameterization allows SRM with stochastic threshold
neurons to learn distinct temporal receptive fields [4, Fig. 5].

###  A.3 Experimental details

Table 2: Hyperparameters used in the different experiments. Parameter \ Experiment | Predictive coding | MNIST | MNIST-DVS
---|---|---|---
$\beta$ | $1$ | $1\text{e}-3$ | $1\text{e}-3$
$\eta$ | $1\text{e}-2$ | $1\text{e}-5$ | $1\text{e}-4$
$N_{X}$ | 20 | 784 | 676
$N_{Y}$ | 10 | 256 | 256
$N_{R}$ | 210 | 784 | 784
Decoder | Linear | MLP/Conv | Conv
Encoder architecture | $N_{Y}$ | $600-N_{Y}$ | $800-N_{Y}$
$\tau_{d}$ | $5$ | $T$ | $T$
$\tau_{e}$ | $5$ | $T$ | $T$
$p$ | $0.2$ | $0.3$ | $0.3$

#### Decoder architectures.

As explained in Sec. 6, we use three types of decoders: a fully-connected
linear model, an MLP, and a convolutional ANN. The MLP comprises two layers,
with architecture $(N_{Y}\times\tau_{d})/2-N_{R}$. The convolutional ANN has
two convolutional layers with ReLU activation and a fully-connected layer,
with architecture
$\frac{T}{2}\mathtt{c}3\mathtt{p}1\mathtt{s}2-20\mathtt{c}3\mathtt{p}1\mathtt{s}2-N_{R}$,
where $X\mathtt{c}Y\mathtt{p}Z\mathtt{s}S$ represents $X$ 1D convolution
filters ($Y\times Y$) with padding $Z$ and stride $S$.

#### Training hidden neurons with SG.

Part of the code used for feedback alignment has been adapted from
https://github.com/ChFrenkel/DirectRandomTargetProjection/, distributed under
the Apache v2.0 license.

Figure 4: MSE for efficient filtering ($\delta=3$) as a function of the window
sizes $\tau_{e}$ and $\tau_{d}$.

#### Efficient predictive coding.

In this experiment, we consider a logistic regression model for the decoder,
i.e.,
$q_{{\bm{w}}^{d}}^{(\tau_{d})}({\bm{r}}_{t}|{\bm{y}}_{\tau_{d}}^{t})=\mathcal{S}({\bm{r}}_{t}|f_{{\bm{w}}^{d}}({\bm{y}}_{\tau_{d}}^{t}))$
where $\mathcal{S}(\cdot)$ is the softmax function and
$f_{{\bm{w}}^{d}}({\bm{x}})={\bm{w}}^{d}{\bm{x}}+{\bm{b}}$ with learnable bias
${\bm{b}}$. We choose $N_{X}=N_{Y}=20$, which gives
$N_{R}=\binom{20}{1}+\binom{20}{2}=210$. Other hyperparameters are precised in
Table 2. Results are averaged over five trials with random initializations of
the networks.

In Fig. 4, we analyze the effect of the window lengths $\tau_{e}$ and
$\tau_{d}$ for filtering (lag $\delta=-3$). The MSE is seen to decrease as the
window sizes $\tau_{e}$ and $\tau_{d}$ grow larger, which allows the SNN and
ANN to access more information when estimating an input sample. The encoding
window is seen to have a prominent effect, with a large improvement in MSE for
$\tau_{e}>3$, i.e., when the system is fed the signal to reconstruct.

Figure 5: Top: Reconstructed MNIST digits. Bottom: Corresponding original
target images.

#### MNIST.

Table 3: Comparison of classification accuracy with various encoding and decoding strategies for naturalization of MNIST digits, with MLP and convolutional ANN decoders. Encoding \ Decoding | Rate (MLP) | VDIB (MLP) | VDIB (Conv)
---|---|---|---
Rate | $60.50\pm 0.32\%$ | $\bm{91.32\pm 0.27\%}$ | $88.73\pm 1.06\%$
Time-to-First-Spike | $86.43\pm 0.84\%$ | $\bm{92.82\pm 0.25\%}$ | $90.72\pm 0.39\%$

The MNIST dataset is freely available online and has been obtained from
https://pytorch.org/vision/stable/_modules/torchvision/datasets/mnist.html. It
consists of $60,000$ training examples from ten classes, and $10,000$ test
examples. Results are averaged over three trials with random initializations
of the networks. In Fig. 5, we show examples of reconstructed digits with time
encoding, and an MLP decoder. This validates that, as well as being able to be
classified with high accuracy, the reconstructions look accurate. We also
provide additional results on accuracy, using a convolutional decoder. This
choice causes a small degradation in the performance of the system, but
alleviates requirements on the GPU memory. Code and weights for the LeNet
network have been obtained from https://github.com/csinva/gan-vae-pretrained-
pytorch and are freely available.

#### MNIST-DVS.

Table 4: Comparison of the reconstruction MSE of MNIST-DVS digits for a convolutional ANN decoder under two decoding schemes (for the output of the SNN). Decoding | Rate | VDIB
---|---|---
MSE | $0.19\pm 0.02$ | $\bm{0.02\pm 0.0006}$

The MNIST-DVS dataset is freely available online from http://www2.imse-
cnm.csic.es/caviar/MNISTDVS.html. Examples have been obtained by displaying
MNIST digits on a screen to a DVS camera. More information on the dataset can
be found from [29].

In Table 4, we provide a quantitative measure of the results obtained for the
reconstruction of MNIST-DVS digits. Specifically, we report the MSE between
reconstructed and target MNIST digits with two decoding strategies. VDIB is
shown to clearly outperform rate decoding.


G_{k}}\left(n_{\max}\cdot\pi_{i}^{\ast}-\hat{w}_{i}\right)_{+}$
$\displaystyle\leavevmode\nobreak\ -{r_{k}}\cdot\sum_{i\in
G_{k}}\left(\frac{(\sum_{i\notin
G_{k}}\pi_{i}^{\ast})\hat{w}_{i}}{(\sum_{i\notin
G_{k}}\widetilde{w}_{i})\pi^{*}_{i}}-1\right)_{+}\cdot\mathbbm{1}\left(\hat{w}_{i}\geq
n_{\max}\cdot\pi_{i}\right)$ $\displaystyle>$
$\displaystyle\leavevmode\nobreak\ \bar{u}_{k}(\widetilde{w}_{G_{k}}).$

This contradicts with the fact that $\widetilde{w}_{G_{k}}$ is the maximizer
of $\bar{u}_{k}$.

#### Step 2.

We show that there exists $\lambda\geq n_{\max}$ such that
$\widetilde{w}=\lambda\cdot\pi^{\ast}$.

Suppose $\exists\leavevmode\nobreak\ k\in[K]$ and $i\in G_{k}$ such that
$\frac{(\sum_{j\notin G_{k}}\pi_{i}^{\ast})\widetilde{w}_{i}}{\sum_{j\notin
G_{k}}\widetilde{w}_{i}}>\pi_{i}^{\ast}$, define
$\epsilon_{i}=\frac{(\sum_{j\notin
G_{k}}\pi_{j}^{\ast})\widetilde{w}_{i}}{(\sum_{j\notin
G_{k}}\widetilde{w}_{j})\pi^{*}_{i}}-1>0$. Let
$\hat{w}_{j}=\begin{cases}\widetilde{w}_{j},\leavevmode\nobreak\
&\text{if}\leavevmode\nobreak\ j\neq i\\\ \frac{(\sum_{j\notin
G_{k}}\widetilde{w}_{j})\pi^{*}_{i}}{\sum_{j\notin
G_{k}}\pi_{j}^{\ast}},\leavevmode\nobreak\ &\text{if}\leavevmode\nobreak\
j=i\end{cases}$. We have

$\displaystyle\leavevmode\nobreak\ \bar{u}_{k}(\widetilde{w}_{G_{k}})$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
-\log\det\left((A^{(k)})^{\top}\left((1+\epsilon_{i})\hat{w}_{i}x_{i}x_{i}^{\top}+\sum_{j\neq
i}\widetilde{w}_{j}x_{j}x_{j}^{\top}\right)^{-1}A^{(k)}\right)-c^{(k)}\cdot\left(\sum_{j\in
G_{k}}\widetilde{w}_{j}\right)$ $\displaystyle\leavevmode\nobreak\
-c^{(k)}\cdot\sum_{j\in
G_{k}}\left(n_{\max}\cdot\pi_{j}^{\ast}-\widetilde{w}_{j}\right)_{+}-{r_{k}}\cdot\sum_{j\in
G_{k}}\left(\frac{(\sum_{l\notin
G_{k}}\pi_{l}^{\ast})\widetilde{w}_{j}}{(\sum_{l\notin
G_{k}}\widetilde{w}_{l})\pi^{*}_{j}}-1\right)_{+}\cdot\mathbbm{1}\left(\hat{w}_{j}\geq
n_{\max}\cdot\pi_{j}\right)$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
-\log\det\left(\left(1+\epsilon_{i}\right)^{-1}(A^{(k)})^{\top}\left(\sum_{j\in
G_{k}}\hat{w}_{j}x_{j}x_{j}^{\top}\right)^{-1}A^{(k)}\right)-{r_{k}}\cdot\epsilon_{i}-c^{(k)}\cdot\left(\sum_{j\in
G_{k}}\hat{w}_{j}\right)$ $\displaystyle\leavevmode\nobreak\
-c^{(k)}\cdot\sum_{j\in
G_{k}}\left(n_{\max}\cdot\pi_{j}^{\ast}-\hat{w}_{j}\right)_{+}-{r_{k}}\cdot\sum_{j\in
G_{k},j\neq i}\left(\frac{(\sum_{l\notin
G_{k}}\pi_{l}^{\ast})\hat{w}_{j}}{(\sum_{l\notin
G_{k}}\widetilde{w}_{l})\pi^{*}_{j}}-1\right)_{+}\cdot\mathbbm{1}\left(\hat{w}_{j}\geq
n_{\max}\cdot\pi_{j}\right)$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
\bar{u}_{k}(\hat{w}_{G_{k}})+{r_{k}}\cdot\log\left(1+\epsilon_{i}\right)-{r_{k}}\cdot\epsilon_{i}$
$\displaystyle<$ $\displaystyle\leavevmode\nobreak\
\bar{u}_{k}(\hat{w}_{G_{k}}).$

This contradicts with the fact that $\widetilde{w}_{G_{k}}$ is the maximizer
of $\bar{u}_{k}$.

As a result, $\frac{(\sum_{i\notin
G_{k}}\pi_{i}^{\ast})\widetilde{w}_{i}}{\sum_{i\notin
G_{k}}\widetilde{w}_{i}}\geq\pi_{i}^{\ast}$ holds for any $i\in[n]$. Examining
this inequality for all $i\in[n]$ yields that
$\widetilde{w}=\lambda\cdot\pi^{\ast}$ for some $\lambda$ which must be
greater than or equal to $n_{\max}$.

#### Step 3.

We show that $\lambda=n_{\max}$.

From the definition of $n_{\max}$, for any $\lambda>n_{\max}$ there exists
$k\in[K]$ such that

$\displaystyle\left(u^{(k)}\circ\mathcal{M}_{\mathrm{eff}}^{(k)}\right)(\lambda\cdot\pi^{\ast})<$
$\displaystyle\leavevmode\nobreak\ v_{\ast}^{(k)}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
\left(u^{(k)}\circ\mathcal{M}_{\mathrm{eff}}^{(k)}\right)(n_{\max}\cdot\pi^{\ast})$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\left(u^{(k)}\circ\mathcal{M}_{\mathrm{eff}}^{(k)}\right)\left((n_{\max}\cdot\pi^{\ast}_{G_{k}},\lambda\cdot\pi^{\ast}_{G_{k}^{c}})\right).$

This means that $\lambda\cdot\pi^{\ast}_{G_{k}}$ is not a NE, yielding a
contradiction. It follows that
$\left(n_{\max}\cdot\pi^{\ast}_{G_{k}}\right)_{k\in[K]}$ is the unique pure
NE.

∎

##  Appendix F Supporting Lemmas

######  Lemma F.1.

We abbreviate $A^{(k)}$ and $A^{(k)}(A^{(k)})^{\top}$ as $A$ and $P$
respectively. Suppose $\mathbf{M}(w)\succ 0$. For any $l\in G_{k}$, we have

$\displaystyle\frac{\partial-\log\det\left(A^{\top}\mathbf{M}(w)^{-1}A\right)}{\partial
w_{l}}=\left\langle x_{l}x_{l}^{\top},\mathbf{M}(w)^{-1}\right\rangle>0.$

Furthermore, fixing $w_{G_{k}^{c}}$, the function
$f(w)=-\log\det\left(A^{\top}\mathbf{M}(w)^{-1}A\right)$ is concave in
$w_{G_{k}}$.

###### Proof.

Indeed, for $l\in G_{k}$

$\displaystyle\leavevmode\nobreak\
\frac{\partial-\log\det\left(A^{\top}\mathbf{M}(w)^{-1}A\right)}{\partial
w_{l}}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
-\left\langle\left(A^{\top}\mathbf{M}(w)^{-1}A\right)^{-1},\frac{\partial
A^{\top}\mathbf{M}(w)^{-1}A}{\partial w_{l}}\right\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
-\left\langle\left(A^{\top}\mathbf{M}(w)^{-1}A\right)^{-1},A^{\top}\frac{\partial\mathbf{M}(w)^{-1}}{\partial
w_{l}}A\right\rangle$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
-\left\langle\left(A^{\top}\mathbf{M}(w)^{-1}A\right)^{-1},A^{\top}\left(\sum_{i,j\in[d]}-\mathbf{M}(w)^{-1}e_{i}e_{j}^{\top}\mathbf{M}(w)^{-1}(x_{l}x_{l}^{\top})_{i,j}\right)A\right\rangle$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
-\left\langle\left(A^{\top}\mathbf{M}(w)^{-1}A\right)^{-1},A^{\top}\left(-\mathbf{M}(w)^{-1}x_{l}x_{l}^{\top}\mathbf{M}(w)^{-1}\right)A\right\rangle$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\left\langle\left(A^{\top}\mathbf{M}(w)^{-1}A\right)^{-1},A^{\top}\mathbf{M}(w)^{-1}AA^{\top}x_{l}x_{l}^{\top}\mathbf{M}(w)^{-1}A\right\rangle$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \left\langle
I,A^{\top}x_{l}x_{l}^{\top}\mathbf{M}(w)^{-1}A\right\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \left\langle
x_{l}x_{l}^{\top}AA^{\top},\mathbf{M}(w)^{-1}\right\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \left\langle
x_{l}x_{l}^{\top},\mathbf{M}(w)^{-1}\right\rangle,$

where the first step uses the fact that $\frac{\partial}{\partial
Y_{i,j}}\log\det Y=[Y^{-1}]_{j,i}$; the third step uses the fact that
$\frac{\partial}{\partial Y_{i,j}}Y^{-1}=-Y^{-1}e_{i}e_{j}^{\top}Y^{-1}$; the
fourth step comes from

$\displaystyle\sum_{i,j\in[d]}e_{i}e_{j}^{\top}(x_{l}x_{l}^{\top})_{i,j}=\sum_{i,j\in[d]}e_{i}e_{i}^{\top}(x_{l}x_{l}^{\top})e_{j}e_{j}^{\top}=x_{l}x_{l}^{\top};$

the fifth and final step use the fact that
$AA^{\top}x_{l}x_{l}^{\top}=x_{l}x_{l}^{\top}$ since $AA^{\top}=P$ is the
projection matrix on $\\{x_{i}\\}_{i\in G_{k}}$.

This establishes the first statement. To show concavity, notice that

$\displaystyle\left[-\mathrm{Hess}_{f}(w)\right]_{i,j}=$
$\displaystyle\leavevmode\nobreak\ -\frac{\partial\left\langle
x_{i}x_{i}^{\top},\mathbf{M}(w)^{-1}\right\rangle}{\partial w_{j}}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \left\langle
x_{i}x_{i}^{\top},\mathbf{M}(w)^{-1}x_{j}x_{j}^{\top}\mathbf{M}(w)^{-1}\right\rangle$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\left(x_{i}^{\top}\mathbf{M}(w)^{-1}x_{j}\right)^{2},\leavevmode\nobreak\
\forall i,j\in G_{k}.$

Therefore $-\mathrm{Hess}_{f}(w)$ is Hadamard product of the positive semi-
definite matrix $\left(x_{i}^{\top}\mathbf{M}(w)^{-1}x_{j}\right)_{i,j\in
G_{k}}$ and itself. It follows from Schur product theorem that the negative
Hessian matrix $-\mathrm{Hess}_{f}$ is symmetric positive semidefinite
everywhere in the domain, which establishes the concavity. ∎

######  Lemma F.2.

For any $w$ such that $\sum_{i=1}^{n}w_{i}x_{i}x_{i}^{\top}$ is non-singular,

$\displaystyle\log\det\left((A^{(k)})^{\top}\left(\sum_{i\in
G_{k}}w_{i}x_{i}x_{i}^{\top}\right)A^{(k)}\right)\leq$
$\displaystyle\leavevmode\nobreak\
\log\det\left((A^{(k)})^{\top}\left(\sum_{i=1}^{n}w_{i}x_{i}x_{i}^{\top}\right)^{-1}A^{(k)}\right)^{-1}.$

Further, if equality holds, then for any $w_{G_{k}}^{\prime}$ we have

$\displaystyle\left(A^{(k)})^{\top}(\sum_{i\in
G_{k}}w_{i}^{\prime}x_{i}x_{i}^{\top})A^{(k)}\right)^{-1}=(A^{(k)})^{\top}\left(\sum_{i\in
G_{k}}w_{i}^{\prime}x_{i}x_{i}^{\top}+\sum_{i\notin
G_{k}}w_{i}x_{i}x_{i}^{\top}\right)^{-1}A^{(k)}.$

###### Proof.

Fix $k\in[K]$. We use shorthand notation $M_{0}=\sum_{i\in
G_{k}}w_{i}x_{i}x_{i}^{\top},M=\sum_{i=1}^{n}w_{i}x_{i}x_{i}^{\top},A=A^{(k)}$.
Let $A=U\Lambda V^{\top}$ denote the Singular Value Decomposition (SVD)
decomposition of $A$, where $U\in\mathbb{R}^{d\times d}$ and
$V\in\mathbb{R}^{r_{k}\times r_{k}}$ are real orthogonal matrices. Then
$M_{0}$ can be written as
$M_{0}=U\begin{pmatrix}D&\mathbf{0}_{r_{k}\times(d-r_{k})}\\\
\mathbf{0}_{(d-r_{k})\times
r_{k}}&\mathbf{0}_{(d-r_{k})\times(d-r_{k})}\end{pmatrix}U^{\top}$ where
$D\in\mathbb{S}^{r_{k}}_{+}$. We write
$U(M-M_{0})U^{\top}=\begin{pmatrix}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{pmatrix}$ where $A_{11}\in\mathbb{R}^{r_{k}\times r_{k}}$.

Notice that

$\displaystyle A^{\top}M_{0}A=$ $\displaystyle\leavevmode\nobreak\
V\Lambda^{\top}U^{\top}U\begin{pmatrix}D&\mathbf{0}_{r_{k}\times(d-r_{k})}\\\
\mathbf{0}_{(d-r_{k})\times
r_{k}}&\mathbf{0}_{(d-r_{k})\times(d-r_{k})}\end{pmatrix}U^{\top}U\Lambda
V^{\top}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
V\begin{pmatrix}\mathbf{I}_{r_{k}\times
r_{k}}&\mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}\begin{pmatrix}D&\mathbf{0}_{r_{k}\times(d-r_{k})}\\\
\mathbf{0}_{(d-r_{k})\times
r_{k}}&\mathbf{0}_{(d-r_{k})\times(d-r_{k})}\end{pmatrix}\begin{pmatrix}\mathbf{I}_{r_{k}\times
r_{k}}\\\ \mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}V^{\top}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ VDV^{\top}.$

We assert that $A^{\top}M^{-1}A\preceq\left(A^{\top}M_{0}A\right)^{-1}$.
Indeed, we have

$\displaystyle\leavevmode\nobreak\
\left(A^{\top}M_{0}A\right)^{1/2}A^{\top}M^{-1}A\left(A^{\top}M_{0}A\right)^{1/2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
VD^{1/2}V^{\top}V\Lambda^{\top}U^{\top}U\begin{pmatrix}D+A_{11}&A_{12}\\\
A_{21}&A_{22}\end{pmatrix}^{-1}U^{\top}U\Lambda V^{\top}VD^{1/2}V^{\top}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
VD^{1/2}\begin{pmatrix}\mathbf{I}_{r_{k}\times
r_{k}}&\mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}\begin{pmatrix}(D+A_{11}-A_{12}A_{22}^{-1}A_{21})^{-1}&\star\\\
\star&\star\end{pmatrix}\begin{pmatrix}\mathbf{I}_{r_{k}\times r_{k}}\\\
\mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}D^{1/2}V^{\top}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
VD^{1/2}\begin{pmatrix}D+A_{11}-A_{12}A_{22}^{-1}A_{21}\end{pmatrix}^{-1}D^{1/2}V^{\top}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ VV^{\top}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \mathbf{I}_{r_{k}\times
r_{k}}$

where the inequality is due to $A_{11}-A_{12}A_{22}^{-1}A_{21}\succeq 0$ since
this is the Schur complement of $M-M_{0}\succeq 0$. As a result,
$A^{\top}M^{-1}A\preceq\left(A^{\top}M_{0}A\right)^{-1}$.

Applying the monotonicity of $\log\det(\cdot)$, we have

$\displaystyle\log\det\left((A^{(k)})^{\top}\left(\sum_{i=1}^{n}w_{i}x_{i}x_{i}^{\top}\right)^{-1}A^{(k)}\right)\leq\log\det\left((A^{(k)})^{\top}\left(\sum_{i\in
G_{k}}w_{i}x_{i}x_{i}^{\top}\right)A^{(k)}\right)^{-1}.$

This establishes the inequality. Further, if equality holds then
$A_{11}-A_{12}A_{22}^{-1}A_{21}=0$. As a result,

$\displaystyle\leavevmode\nobreak\ A^{\top}M^{-1}A$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
V\Lambda^{\top}U^{\top}U\begin{pmatrix}D+A_{11}&A_{12}\\\
A_{21}&A_{22}\end{pmatrix}^{-1}U^{\top}U\Lambda V^{\top}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ V\begin{pmatrix}\mathbf{I}_{r_{k}\times
r_{k}}&\mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}\begin{pmatrix}(D+A_{11}-A_{12}A_{22}^{-1}A_{21})^{-1}&\star\\\
\star&\star\end{pmatrix}\begin{pmatrix}\mathbf{I}_{r_{k}\times r_{k}}\\\
\mathbf{0}_{r_{k}\times(d-r_{k})}\end{pmatrix}V^{\top}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ VD^{-1}V^{\top}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \left(A^{\top}M_{0}A\right)^{-1}.$

Notice that the validity above argument does not depend on $w_{G_{k}}$ and in
fact holds for any $w_{G_{k}}^{\prime}$. This completes the proof.

∎

######  Claim F.3 (Kakutani’s fixed point theorem).

Consider a set-valued function $F:\mathcal{C}\rightarrow 2^{\mathcal{C}}$ over
convex and compact set $\mathcal{C}$ such that (i) $F$ has closed graph; (ii)
$F(x)$ is non-empty and convex for any $x\in\mathcal{C}$, then there exists a
fixed point $x$ such that $x\in F(x)$.

######  Claim F.4 (Monotonicity of determinant).

Suppose $A$ and $B$ are two symmetric matrices such that $A\succeq B\succ 0$,
then $\det A\geq\det B$.

######  Claim F.5 (Concavity of log-determinant function).

Suppose $A$ and $B$ are two symmetric positive semidefinite matrices such that
$A\succeq B\succ 0$, then $\log\det(\lambda A+(1-\lambda)B)\geq\lambda\log\det
B+(1-\lambda)\log\det A$ holds for any $\lambda\in(0,1)$.

The following important result of Kiefer and Wolfowitz established the
equivalence of the D-optimal design and G-optimal design.

######  Theorem F.6 (General equivalence theorem of G-optimal design [33]).

Assume $\mathrm{span}(\mathcal{X})=\mathbb{R}^{d}$. The followings are
equivalent:

  * •

$\pi^{\ast}=\underset{\pi\in\Delta(\mathcal{X})}{\arg\max}\leavevmode\nobreak\
\log\det\mathbf{M}(\pi)$;

  * •

$\pi^{\ast}=\underset{\pi\in\Delta(\mathcal{X})}{\arg\min}\leavevmode\nobreak\
\max_{x\in\mathcal{X}}x^{\top}\mathbf{M}(\pi)^{-1}x$;

  * •

$\max_{x\in\mathcal{X}}x^{\top}\mathbf{M}(\pi)^{-1}x=d$.