Signal Coding and Perfect Reconstruction using Spike Trains

TL;DR

This paper introduces a new mathematical framework for encoding and perfect reconstruction of signals using biologically plausible spike trains, leveraging convex optimization and stochastic gradient descent.

Contribution

It presents a novel model for signal coding and reconstruction with formal conditions for perfect recovery, applicable to neural spike train analysis.

Findings

Framework enables perfect signal reconstruction under certain conditions

Convex optimization effectively reconstructs signals from spike trains

Simulation results demonstrate the method's robustness and efficacy

Abstract

In many animal sensory pathways, the transformation from external stimuli to spike trains is essentially deterministic. In this context, a new mathematical framework for coding and reconstruction, based on a biologically plausible model of the spiking neuron, is presented. The framework considers encoding of a signal through spike trains generated by an ensemble of neurons via a standard convolve-then-threshold mechanism. Neurons are distinguished by their convolution kernels and threshold values. Reconstruction is posited as a convex optimization minimizing energy. Formal conditions under which perfect reconstruction of the signal from the spike trains is possible are then identified in this setup. Finally, a stochastic gradient descent mechanism is proposed to achieve these conditions. Simulation experiments are presented to demonstrate the strength and efficacy of the framework

Equations

$$\int X(\tau )K^{j_i}(t_i-\tau )d\tau =T^{j_i}(t_i)$$

$$T^{j}(t)=\begin{cases*}C^{j},&$t-\delta_{j}>t_{l}^{j}(t)$\\ M^{j}-\frac{(t-t_{l}^{j}(t))(M^{j}-C^{j})}{\delta_{j}},&$t-\delta_{j}\leq t_{l}^{j}(t)$\\ \end{cases*}$$

$$X^{\ast }(t)=\underset{\tilde{X}}{\mathrm{argmin}}||\tilde{X}(t)|{|}_2^2$$

$$\text{s.t.}\int \tilde{X}(\tau )K^{j_i}(t_i-\tau )d\tau =T^{j_i}(t_i);1\le i\le N$$

$$\mathcal{G}=\{X(t)|X(t)\in \mathcal{F},X(t)=\sum _{p=1}^N{\alpha }_p{K}^{j_p}(t_p-t),j_p\in \{1,...,n\},{\alpha }_p\in R,t_p\in R^{+},N\in Z^{+}\}$$

$$X^{\ast }(t)=\sum _{i=1}^N{\alpha }_i{K}^{j_i}(t_i-t)$$

$${\forall }_{1\le i\le N};\int \sum _{k=1}^N{\alpha }_k{K}^{j_k}(t_k-t)K^{j_i}(t_i-t)d\tau =T^{j_i}(t_i)$$

$${\forall }_{1\le i\le N};\sum _{k=1}^N{P}_{ik}{\alpha }_k=b_i$$

$$P\alpha =T$$

$$\tilde{X}(t)=\sum _{i=1}^N{a}_i{K}^{j_i}(t_i-t)$$

$$||X(t)-X^{\ast }(t)||\le ||X(t)-\tilde{X}(t)||$$

$$||X(t)-\tilde{X}(t)||=||\underset{\text{A}}{\underbrace{X(t)-X^{\ast }(t)}}+\underset{\text{B}}{\underbrace{X^{\ast }(t)-\tilde{X}(t)}}||$$

$$\text{First, }⟨A,K^{j_i}(t_i-t)⟩=⟨X(t),K^{j_i}(t_i-t)⟩-⟨X^{\ast }(t),K^{j_i}(t_i-t)⟩,\forall i\in \{1,2,..,N\}$$

$$=T^{j_i}(t_i)-T^{j_i}(t_i)=0\text{ (Using the constraints in (refoptimizationproblem) \& (refspikeConstraint))}$$

$$\text{Second, }⟨A,B⟩=⟨A,\sum _{i=1}^N({\alpha }_i-a_i)K^{j_i}(t_i-t)⟩\text{ (By Lemma 1 X∗(t)=∑i=1NαiKji(ti−t))}$$

$$=\sum _{i=1}^N({\alpha }_i-a_i)⟨A,K^{j_i}(t_i-t)⟩=0$$

$$⟹||X(t)-\tilde{X}(t)|{|}^2=||A+B|{|}^2=||A|{|}^2+2⟨A,B⟩+||B|{|}^2$$

$$=||A|{|}^2+||B|{|}^2\ge ||A|{|}^2=||X(t)-X^{\ast }(t)|{|}^2$$

$$⟹||X(t)-\tilde{X}(t)||\ge ||X(t)-X^{\ast }(t)||\square$$

$$X(t)=\sum _{p=1}^N{\alpha }_p{K}^{j_p}(t_p-t)({\alpha }_p\in R,t_p\in R^{+},N\in Z^{+})$$

$$⟨X(t),K^{j_p}(t_p-t)⟩=T^{j_p}(t_p)\forall p\in \{1,...,N\}$$

$$||X(t)-X^{\ast }(t)|{|}_2\le ||X(t)-X(t)|{|}_2=0\square$$

$$X(t)=\sum _{i=1}^N{\alpha }_i{f}_{p_i}(t-t_i)$$

$$\text{nicefrac}||X(t)-X_{hyp}(t)|{|}_2^2||X(t)|{|}_2^2\le (\delta +C\Delta )\text{nicefrac}(1+x_{max})(1-\eta )$$

$$X^{\ast }(t)$$

$$K^{j_i}$$

$$X(t)$$

$$\{(t_i^{\prime },K^{j_i})|i\in \{1,...,N\}\}$$

$$X(t)$$

$$M$$

$$X^{\ast }(t)=\sum _{i=1}^N{\alpha }_i{K}^{j_i}(t-t_i^{\prime })+\sum _{k=1}^M\widetilde{{\alpha }_k}{K}^{\widetilde{j_k}}(t-\widetilde{t_k})$$

$${\alpha }_i$$

$$T_i$$

$$<X(t),K^{j_i}(t-t_i^{\prime })>=T_i\forall i\in \{1,2,...,N\}$$

$$X_{hyp}(t)$$

$$X_{hyp}(t)=\sum _{i=1}^N{a}_i{K}^{j_i}(t-t_i^{\prime }),a_i\in R$$

$$s.t.<X_{hyp}(t),K^{j_i}(t-t_i^{\prime })>=T_i\forall i$$

$$a_i$$

$$t_i^{\prime }$$

$$X_{hyp}(t)$$

$$X(t)$$

Full text

1 Introduction

Spike based encoding of sensory stimuli is a hallmark of biological systems. It is now well-established that the coding of continuous time sensory signals in spike trains is a complex and diverse phenomenon, and is fairly deterministic in many animal sensory pathways[12, 5, 14, 19, 9, 13]. Spike train representations, when sparse, are not only intrinsically energy efficient, but can also facilitate computation at later stages of processing[7, 8]. In their seminal work, Olshausen and Field [15] showed how efficient codes can arise from learning sparse representations of natural stimulus statistics, resulting in striking similarities with observed biological receptive fields. Smith and Lewicki [10, 3] likewise showed that auditory filters could be estimated by training a population spike code model with natural sounds. These studies, by and large, fall under the general framework of dictionary learning: identifying an over-complete dictionary $\{\phi_{j}|j=1\ldots m\}$ such that each stimulus $s_{i}$ in an ensemble $\{s_{i}|i=1\ldots n\}$ can be represented as $s_{i}=\sum_{j=1}^{m}\alpha_{j}\phi_{j}$ where the vector of coefficients $\alpha_{j}$ is sparse. The studies only made passing reference to how the $\alpha_{j}$’s may be derived (e.g. matching pursuit [11]) or even be represented (e.g. local population of neurons spiking probabilistically proportional to $\alpha_{j}$ in [3]). Lacking clearly specified plausible neural implementations, the extent to which the proposed schemes underlie biological sensory processing therefore remained unclear. To remedy this, several subsequent learning techniques based on biologically plausible models of spiking neurons have been proposed. For example, [20] developed a biophysically motivated spiking neural network which for the first time predicted the full diversity of V1 simple cell receptive field shapes when trained on natural images. Elsewhere [17] presented a rate encoded spiking neural network of integrate-and-fire neurons demonstrating convergence to nearly optimal encodings.

Although these results signify substantial progress, the classical signal processing question of what class of signals support perfect or approximate reconstruction when coded using spike trains, remains to be fully resolved. Admittedly, the very coarse $\Sigma\Delta$ quantization of bandlimited signals investigated in [6] does amount to a spike train representation. However, due to the classical nature of its framework, not only is biological plausibility not a concern, but also coding is explored in the oversampled regime. Along similar lines, [4] has explored the spike generating mechanism of the neuron as an oversampling, noise shaping analog-to-digital converter.

Here we present a new framework for coding and reconstruction that begins with a biologically plausible coding mechanism which is a superset of the standard leaky integrate-and-fire mechanism. Reconstruction is first formulated as an optimization that minimizes the energy of the reconstructed signal subject to consistency with the spike train, and then solved in closed form. We then identify a general class of signals for which reconstruction is provably perfect under certain conditions. Surprisingly, the result instantiates a version of Barlow’s “efficient coding hypothesis” [1], which posits that the coding strategy of sensory neurons should be adapted to the statistics of the stimuli in an animal’s natural environment. We present a stochastic gradient descent mechanism to achieve these conditions, and close with simulation experiments that demonstrate the efficacy of the framework. The rest of the paper is structured as follows. In Sections 2 and 3 we introduce the coding and decoding frameworks. Section 4 identifies the class of signals for which perfect reconstruction is achievable if certain conditions are met. A learning scheme to achieve these conditions is presented in Section 6, followed by simulation results in Section 7. We conclude in Section 8.

2 Coding

The general class of deterministic transformations (i.e., the set of all nonlinear operators) from continuous time signals to spike trains is difficult to characterize because the space of all spike trains does not lend itself to a natural topology that is universally agreed upon. The result is that simple characterizations, such as the set of all continuous operators, can not be posited in a manner that has general consensus. To resolve this issue, we take a que from biological systems. In most animal sensory pathways, external stimulus passes through a series of transformations before being turned into spike trains[18]. For example, visual signal in the retina is processed by multiple layers of non-spiking horizontal, amacrine and bipolar cells, before being converted into spike trains by the retinal ganglion cells. Accordingly, we can consider the set of transformations that pass via an intermediate continuous time signal which is then transformed into a spike train through a simple stereotyped mapping where spikes mark threshold crossings. The complexity of the operator now lies in the mapping from the continuous time input signal to the continuous time intermediate signal. Since any time invariant, continuous, nonlinear operator with fading memory can be approximated by a finite Volterra series operator[2], this general class of nonlinear operators from continuous time signals to spike trains can be modeled as the composition of a finite Volterra series operator and a neuronal thresholding operation to generate a spike train. In our encoding model the simplest subclass of these transformations is considered: the case where the Volterra series operator has a single causal, bounded-time, linear term, the output of which is composed with a thresholding operation of a potentially time varying threshold. The overall operator from the input signal to the spike train remains nonlinear due to the thresholding operation.

Formally, we assume the input signal $X(t)$ to be a bounded continuous function in the interval $[0,L]$ for some $L\in R^{+}$, i.e., we are interested in the class of input signals $\mathcal{F}=\{X(t)|X(t)\in C[0,L]\}$. Since the framework involves signal snippets of arbitrary length, this choice of $L$ is without loss of generalization. We assume an ensemble of convolution kernels $K=\{K^{j}|j\in Z^{+},j\leq n\}$, consisting of $n$ kernels $K^{j},j=1,\ldots,n$. We assume that $K^{j}(t)$ is a continuous function on a bounded time interval $[0,T]$, i.e. $\forall j\in\{1,\ldots,n\},K^{j}(t)\in C[0,T]$ for some $T\in R^{+}$. Finally, we assume that $K^{j}$ has a time varying threshold denoted by $T^{j}(t)$.

The ensemble of convolution kernels $K$ encodes a given input signal $X(t)$ into a sequence of spikes $\{(t_{i},K^{j_{i}})\}$, where the $i^{th}$ spike is produced by the $j_{i}^{th}$ kernel $K^{j_{i}}$ at time $t_{i}$ if and only if:

[TABLE]

We assume that the time varying threshold $T^{j}(t)$ of the $j{th}$ kernel remains constant at $C^{j}$ until that kernel produces a spike, at which time an after-hyperpolarization potential (ahp) kicks in to raise the threshold to a high value $M^{j}\gg C^{j}$, which then drops back linearly to its original value within a refractory period $\delta_{j}$. Formally, the threshold function $T^{j}(t)$ of the $j^{th}$ kernel is given by:

[TABLE]

Where $t_{l}^{j}(t)$ denotes the time of the last spike generated by $K^{j}$ prior to time $t$.

Notably, apart from the contribution due to the ahp, we have considered the threshold of the $j^{th}$ kernel (alternately called a neuron) to be a constant $C^{j}$ in our model. This is in contrast to real neurons that are known to vary their thresholds through homeostasis.

3 Decoding

The objective of the decoding module is to reconstruct the original signal from the encoded spike trains. Considering the prospect of the invertibility of the coding scheme, we seek a signal that satisfies the same set of constraints as the original signal when generating all spikes apropos the set of kernels in ensemble $K$. Recognizing that such a signal might not be unique, we choose the reconstructed signal as the one with minimum $L2$-norm. Formally, the reconstruction (denoted by $X^{*}(t)$) of the input signal $X(t)$ is formulated to be the solution to the optimization problem:

[TABLE]

where $\{(t_{i},K^{j_{i}})|i\in\{1,...,N\}\}$ is the set of all spikes generated by the encoder.

Why L2 minimization?

The choice of $L2$ minimization as the objective of the reconstruction problem can only be weakly justified at the current juncture. The perfect reconstruction theorem that follows provides the strong justification. As it stands, the $L2$ minimization objective is in congruence with the dictum of energy efficiency in biological systems. The assumption is that, of all signals, the one with the minimum energy that is consistent with the spike trains is desirable. Secondly, an $L2$ minimization in the objective of (3) reduces the convex optimization problem to a solvable linear system of equations as shown in Lemma 1. Later we shall show that $L2$-minimization has the surprising benefit of recovering the original signal perfectly under certain conditions.

4 Signal Class for Perfect Reconstruction

To establish the effectiveness of the described coding-decoding model, we have to evaluate the accuracy of reconstruction over a class of input signals. We observe that in general the encoding of continuous time signals into spike trains is not a one-to-one map; the same set of spikes can be generated by different continuous time signals with appropriate changes in amplitudes so as to result in the same convolved values at the spike times. Naturally, with a finite and fixed ensemble of kernels $K$, one cannot achieve perfect reconstruction for the general class of signals $\mathcal{F}$ as defined in Section 2. We now restrict ourselves to a subset $\mathcal{G}$ of the original class $\mathcal{F}$ as defined below and address the question of reconstruction accuracy.

[TABLE]

Essentially $\mathcal{G}$ consists of all linear combinations of arbitrarily shifted inverted kernel functions. $N$ is bounded above by the total number of spikes that the ensemble $K$ can generate over $[0,L]$. For this restricted class of signals the perfect reconstruction theorem is presented below. The theorem is proved with the help of two lemmas.

Perfect Reconstruction Theorem:

Let $X(t)\in\mathcal{G}$ be an input signal. Then for appropriately chosen time-varying thresholds of the kernels, the reconstruction, $X^{*}(t)$, resulting from the proposed coding-decoding framework is accurate with respect to the $L2$ metric, i.e., $||X^{*}(t)-X(t)||_{2}=0$.

Lemma1:

The solution $X^{*}(t)$ to the reconstruction problem given by (3) can be written as:

[TABLE]

where the coefficients $\alpha_{i}\in R$ can be uniquely solved from a system of linear equations if the shifted kernel functions $K^{j_{i}}(t_{i}-t)$ are linearly independent.

Proof:

Application of the Representer Theorem [16] on (3) directly results in:

[TABLE]

where the $\alpha_{i}$’s are real valued coefficients. In essence, the reconstructed signal $X^{*}(t)$ becomes a summation of the kernels, shifted to their respective times of generation of spikes, scaled by appropriate coefficients. Plugging (6) into the constraints (3) gives:

[TABLE]

Setting $b_{i}=T^{j_{i}}(t_{i})$ and $P_{ik}=\int K^{j_{k}}(t_{k}-\tau)K^{j_{i}}(t_{i}-\tau)d\tau$ results in:

[TABLE]

Equation (7) defines a system of $N$ equations in $N$ unknowns of the form:

[TABLE]

where $\alpha=\langle\alpha_{1},...,\alpha_{N}\rangle^{T}$, $T=\langle T^{j_{1}}(t_{1}),...,T^{j_{N}}(t_{N})\rangle^{T}$ and $P$ is an $N\times N$ matrix with elements $P_{ik}=\int K^{j_{k}}(t_{k}-\tau)K^{j_{i}}(t_{i}-\tau)d\tau$. Clearly $P$ is the Gramian Matrix of the shifted kernels $\{K^{j_{i}}(t_{i}-t)|i\in 1,2,...,N\}$ in the Hilbert space with the standard inner product. Hence $\alpha$ has a unique solution if and only if $P$ is invertible. Finally, the Gram matrix $P$ is invertible if and only if the set of vectors $\{K^{j_{i}}(t_{i}-t)|i\in 1,2,...,N\}$ in Hilbert space are linearly independent. $\hfill\Box$

Import:

The goal of the optimization problem is to find the best object in the feasible set. However, the application of the Representer Theorem, as above, converts the constraints into a determined system of unknowns and equations, turning the feasible set into a single point, effectively changing the optimization problem into a solvable system that results in a closed form solution for the $\alpha_{i}$’s. This implies that instead of solving (3), we can solve for the reconstruction from $X^{*}(t)=\sum_{i=1}^{N}\alpha_{i}K^{j_{i}}(t_{i}-t)$, where $\alpha_{i}$ is the i-th element of $\alpha=P^{-1}T$.

Lemma2:

Let $X^{*}(t)$ be the reconstruction of an input signal $X(t)$ and $\{(t_{i},K^{j_{i}})\}_{i=1}^{N}$ be the set of spikes generated. Then, for any arbitrary signal $\tilde{X}(t)$ within the span of $\{K^{j_{i}}(t_{i}-t)|i\in\{1,2,...,N\}\}$, i.e., the set of shifted inverted kernels at respective spike times, given by

[TABLE]

the following inequality holds:

[TABLE]

Proof:

[TABLE]

Import:

The implication of the above lemma is quite remarkable. The objective defined in (3) chooses a signal with minimum energy satisfying the constraints, deemed the reconstructed signal. However as the lemma demonstrates, this signal also has the minimum * error* with respect to the input signal in the span of the shifted kernels. This signifies that our choice of the objective in the decoding module not only draws from biologically motivated energy optimization principles, but also performs optimally in terms of reconstructing the original input signal within the span of the appropriately shifted spike generating kernels.

Exploring further, for a given input signal $X(t)$ if $S_{1}$ and $S_{2}$ are two sets of spike trains where $S_{1}\subset S_{2}$ produced by two different kernel ensembles, the second a superset of the first, then Lemma 2 further implies that the reconstruction due to $S_{2}$ is at least as good as the reconstruction due to $S_{1}$ because the reconstruction due to $S_{1}$ is in the span of the shifted kernel functions of $S_{2}$ as $S_{1}\subset S_{2}$. This immediately leads to the conclusion that for a given input signal the more kernels we add to the ensemble the better the reconstruction, provided the kernels maintain linear independence.

Proof of the Theorem:

The proof of the theorem follows directly from Lemma 2. Since the input signal $X(t)\in\mathcal{G}$, let $X(t)$ be given by the equation below:

[TABLE]

Assume that the time varying thresholds of the kernels in our kernel ensemble $K$ is set in such a manner that the below conditions are satisfied:

[TABLE]

i.e., each of the kernels $K^{j_{p}}$ at the very least produces a spike at time $t_{p}$ against $X(t)$ (regardless of other spikes at other times). Clearly then $X(t)$ lies in the span of the appropriately shifted and inverted response functions of the spike generating kernels. Applying Lemma 2 it follows that:

[TABLE]

Import:

In addition to demonstrating the potency of the coding-decoding scheme, this theorem frames Barlow’s efficient coding hypothesis [1]—that the coding strategy of sensory neurons be adapted to the statistics of the stimuli—in mathematically concrete terms. Going by the theorem, the spike based encoding necessitates the signals to be in the span of the encoding kernels for perfect reconstruction. Inverting the argument, kernels must learn to adapt to the basis elements that generate the signal corpora for superior reconstruction. A practical challenge on which the reconstruction accuracy depends, as indicated by the theorem, is whether we can generate spikes at the correct temporal locations. One way to tackle this problem, as adopted in our experiments in Section 7, is to set the initial threshold values $C^{j}$s and the refractory periods $\delta_{j}$s to be low in our simple thresholding model (2). This then ensures that the real spikes deviate from their desired locations by at most $\delta_{j}$ time. Since the decoding is a continuous transformation, this changes the reconstruction accuracy by a small amount as is confirmed by our experiments in Section 7.

5 Approximate Reconstruction of a Signal:

The perfect reconstruction theorem 2 essentially stipulates the exact conditions under which exact recovery of a signal is feasible in the proposed framework. But in real applications it could be challenging to meet those conditions for any arbitrary signal. For example, we may not be able to generate the spikes in exact times or the signal may not perfectly fit in the subspace generated from a finite bag of kernels, $\mathcal{G}$ as stated in the theorem 2. The goal of the approximate reconstruction theorem is to give a lower bound on the reconstruction error when the conditions are reasonably relaxed.

Approximate Reconstruction Theorem:

Let $X(t)\in L^{2}([0,L]),\text{for some }L\in R^{+}$, a square integrable continuous time signal in the interval $[0,L]$ be an input to our proposed framework. If the following assumptions are true:

• •

$X(t)$ can be realized as a linear combination of some component signals as below:

[TABLE]

where $p_{i}\in Z^{+}$, and components are chosen from a possible infinite set of functions of unit $L^{2}$ norm: $\{f_{i}(t)|i\in Z^{+},||f_{i}(t)||_{2}=1\}$ i.e. each component function $f_{i}(t)$ is normalized to have unit $L^{2}$ norm, and $\alpha_{i}$ are the real coefficients.

Also there $\exists$ a kernel $K^{j_{i}}$ from our bag of kernels K, such that $f_{p_{i}}$ is very close to $K^{j_{i}}(t)$ in $L^{2}$ norm i.e. $\exists\delta\in R^{+}st||f_{p_{i}}(t)-k^{j_{i}}(t)||_{2}<\delta\forall i\in\{1,...,N\}$.

• •

Also assume that when $X(t)$ is encoded using our framework 2, each one of these fitting kernels $k^{j_{i}}$ produces a spike at time $t_{i}^{\prime}$ such that $|t_{i}-t_{i}^{\prime}|<\Delta\forall i$.

• •

Assume the response function of each kernel in our kernel bag satisfies a Lipschutz type of condition as follows, i.e. $\exists C\in Rs.t.||k^{j}(t)-k^{j}(t-\Delta t)||_{2}\leq C|\Delta t|\forall\Delta t\in R,\forall j$.

• •

Finally, we assume the shifted component functions also satisfy a frame bound type of condition as follows: $\sum_{k\neq i}<f_{p_{i}}(t-t_{i}),f_{p_{k}}(t-t_{k})>\leq\eta\forall i\in{1,...,N}$

Then, the noise to signal ratio in reconstruction of $X(t)$ resulting from the proposed coding-decoding framework is bounded. Specifically, the following inequality is satisfied:

[TABLE]

where $x_{max}$ is a positive number $\in[0,N-1]$ that depends on the maximum overlap of the support of component functions $f_{p_{i}}(t-t_{i})$.

Proof of the Theorem:

[TABLE]

[TABLE]

[TABLE]

6 Kernel Adaptation

As demonstrated by the perfect reconstruction theorem, over a given class of input signals ${F}$, reconstructions are perfect when the kernels from the ensemble $K$ match the (unknown) underlying components that generate the signals in ${F}$. With the goal of improving the quality of the reconstruction, we now propose a gradient descent based method that incrementally changes the kernels to decrease reconstruction error, which then indirectly induces the kernels to fit the unknown components from which signals in ${F}$ are generated.

Informally speaking, there are two effects that a perturbation of a kernel has on the reconstructed signal: (i) perturbation of the kernel $K^{j_{i}}$ directly impacts the reconstruction $\sum\alpha_{i}K^{j_{i}}$, and (ii) perturbation of the kernel perturbs the spike times which then have an impact on the locations at which the kernels are situated to be summed. The perturbation of the spike times also incurs a “domino effect” on future spike times via their ahps.

To be able to apply gradient descent on the kernels, we consider the kernels in their parametric forms, i.e., the kernel response function of $K^{j}$ is assumed to be of the form, $K^{j}(t)=K^{j}(t;\{\theta^{j}_{c}\})$, where $\{\theta^{j}_{c}\}$ are the free parameters that govern the shape of the kernel $K^{j}$. In this context, our goal is to derive an update rule for the free parameters $\theta^{j}_{c}$s.

Formally, let us consider a family of signals, $F$, and let $X^{*}(t)$ be the reconstruction of $X(t)$ for $X(t)\in F$. The goal is to arrive at the optimal set of kernels that minimize the expected reconstruction error $\int||X(t)-X^{*}(t)||^{2}dF_{X}$ using gradient descent, where $F_{X}$ is the cumulative probability distribution over $F$ at $X(t)$. Shifting to a stochastic gradient descent framework, for a randomly selected signal $X(t)\in F$, the reconstruction error of $X(t)$, using (6), is defined as:

[TABLE]

Derivative of $t_{i}$.

First, we notice that each spike is generated by a certain neuron, and only the $\theta^{j}_{c}$s of the kernel response of that neuron have an effect on the time of the spike. In other words, if $j_{i}\neq k$ then $\forall{c};\nicefrac{{\partial t_{i}}}{{\partial\theta_{c}^{k}}}=0$. For the other case, we can calculate $\nicefrac{{\partial t_{i}}}{{\partial\theta_{c}^{j_{i}}}}$ as follows: The initial condition for the generation of a spike at time $t_{i}$ by the $j_{i}$th kernel is given by equation(1). Now consider that out of all the free parameters, we only perturb $\theta_{c}^{j_{i}}$ to $\theta_{c}^{j_{i}}+\Delta\theta_{c}^{j_{i}}$. This will in turn perturb all the spike times for the $j_{i}$th kernel. Let spike time $t_{i}$ be shifted to $t_{i}+\Delta t_{i}$. Let us also assume that $t_{l}$ is the last spike produced by $K^{j_{i}}$ prior to $t_{i}$, and $t_{l}$ is correspondingly shifted to $t_{l}+\Delta t_{l}$. Under this perturbed scenario, we can rewrite (1) as:

[TABLE]

Combining (1) and (28), using Taylor series approximation, ignoring all second and higher order terms, using the form of the simple threshold function $T^{j_{i}}(t)$ as in (2) and setting $\lim{\Delta\theta_{c}^{j_{i}}\to 0}$, the derivative of $t_{i}$ with respect to a particular kernel parameter $\theta_{c}^{j_{i}}$ is obtained as:

[TABLE]

It is important to note that the formula for $\frac{\partial t_{i}}{\partial\theta_{c}^{j_{i}}}$ involves derivative of previous spike times, i.e., $\frac{\partial t_{l}}{\partial\theta_{c}^{j_{i}}}$s and hence can be computed recursively. This is the “domino effect” referred to earlier.

Derivative of $E_{X}$.

Using equation (27), we can find the derivative of the error functional with respect to the free parameters as:

[TABLE]

[TABLE]

Using the constraint in equation (3), we see that the first term above goes to [math], leaving the other two remaining terms. We emphasize that this corresponds to the fact that $\frac{\partial E_{X}}{\partial\theta_{c}^{j}}$ depends upon the $\alpha_{i}$s but not on the $\frac{\partial\alpha_{i}}{\partial\theta_{c}^{j}}$.

For a chosen set of kernels, given in a particular parametric form, equations (30) and (29) hold all the information necessary to apply stochastic gradient descent on the $\theta_{c}^{j}$s to minimize the expected reconstruction error.

7 Experiments

The experiments were targeted toward establishing the effectiveness of the proposed learning technique. To achieve this, we first chose an ensemble of $n$ kernels $K$ in a specific parametric form. Here, we used B-splines of order 3 to construct our kernels. Specifically, for kernel $K^{j}\in K$, its response function was of the form:

[TABLE]

where, $B(t)$ is the standard continuously differentiable B-spline function of order 3. In our experiments $\beta^{j}_{c}$s $\in R$ were chosen as the free parameters for the kernels and were randomly initialized. The rest of the parameters were kept fixed throughout the experiments. Next, to construct a class of input signals ${F}$ another ensemble of $n$ functions ${J}$ built out of similar B-splines but with different values of the free parameters was chosen. Specifically an $f^{j}(t)\in{J}$ was given by:

[TABLE]

Using this ensemble of functions $J$, the class of input signals ${F}$ were constructed as a linear combinations of randomly shifted $f^{j}(t)$s, i.e., any $X(t)\in{F}$ was of the form:

[TABLE]

In the experiments, at each step a random sample $X(t)$ from the class ${F}$ was chosen. We then made a gradient update to the free parameters $\beta^{j}_{c}$s of the kernels $K^{j}$s as described in Section 6. If the proposed learning were effective, the reconstruction error would drop rapidly resulting in near perfect reconstruction. We found this to be the case in all our experiments, as reported below.

In repeated experiments, and across varied number ($n$) of differently initialized kernels, we saw steady drop in reconstruction error and with sufficiently large number of learning iterations reconstruction became perfect. Here, we have reported results of sample experiments with 1,5 and 10 kernels in Figure  1. In each case at least $100K$ iterations of learning were executed and after each $1000$ steps of training, kernels were extracted to be run against $1000$ randomly generated new input samples for reporting testing error. The mean reconstruction error on the test set and its standard deviation have been reported at regular intervals using error-bars in the same Figure  1.

One of the primary reasons behind choosing B-splines was its universal approximation property in that any arbitrary bounded-time continuous function can be approximated by B-splines. In that sense a demonstration of learning with B-spline kernels establishes generality. For learning to work we had to tune the parameters properly (e.g. setting ahp refractory period, initial thresholds to low values, putting a reasonable learning rate, etc.) all of which can be found in configurations in our java-based implementation. All experiments were run on standard 8GB machines with quad-core processors taking reasonable time($\sim 1$sec per iteration).

8 Conclusion

In this paper, we have formulated a framework that identifies the precise conditions under which continuous time signals can be represented using an ensemble of spike trains, from which the signal can be recovered perfectly. Although aligned in their goals, this framework is very different from that investigated in Niquist-Shanon theory. The primary difference between the two lies in their respective modes of representation/coding. Instead of sampling the value of a function at uniform or non-uniform prespecified sample points the new coding scheme reports the (non-uniform) sample points where the function takes specific convolved values. A natural extension of this framework addresses approximate reconstruction, both when there is noise injected in the spike timings as well as in the kernel functions. Coding is intimately related to compression and our experimental results indicate great potential in this regard. The simulation source code is available at: bitbucket.org/crystalonix/oldsensorycoding.git.