Sparse Functional Identification of Complex Cells from Spike Times and the Decoding of Visual Stimuli

TL;DR

This paper presents a novel sparse functional identification method for complex cells that efficiently decodes visual stimuli from spike data, outperforming existing models and algorithms.

Contribution

It introduces a rank minimization framework for decoding and identifying complex cell responses, establishing a duality that enhances evaluation and performance.

Findings

Significantly reduces sampling measurements needed for decoding

Outperforms generalized quadratic and spike-triggered covariance models

Provides efficient algorithms for low-rank dendritic stimulus identification

Abstract

We investigate the sparse functional identification of complex cells and the decoding of visual stimuli encoded by an ensemble of complex cells. The reconstruction algorithm of both temporal and spatio-temporal stimuli is formulated as a rank minimization problem that significantly reduces the number of sampling measurements (spikes) required for decoding. We also establish the duality between sparse decoding and functional identification, and provide algorithms for identification of low-rank dendritic stimulus processors. The duality enables us to efficiently evaluate our functional identification algorithms by reconstructing novel stimuli in the input space. Finally, we demonstrate that our identification algorithms substantially outperform the generalized quadratic model, the non-linear input model and the widely used spike-triggered covariance algorithm.

Equations

$$u_1(t)=\sum _{l_t=-L_t}^{L_t}{c}_{l_t}{e}_{l_t}(t),$$

$$e_{l_t}(t)=\frac{1}{\sqrt{S_t}}\exp \left(\frac{j{l}_t{\Omega }_t}{L_t}t\right).$$

$$K_1(t;t^{\prime })=\sum _{l_t=-L_t}^{L_t}{e}_{l_t}(t-t^{\prime }).$$

$$u_2(t_1;t_2)=\sum _{l_{t_1}=-L_t}^{L_t}\sum _{l_{t_2}=-L_t}^{L_t}{d}_{l_{t_1}{l}_{t_2}}{e}_{l_{t_1}}(t_1)\cdot e_{l_{t_2}}(t_2)$$

$$K_2(t_1,t_2;t_1^{\prime },t_2^{\prime })=\sum _{l_{t_1}=-L_t}^{L_t}\sum _{l_{t_2}=-L_t}^{L_t}{e}_{l_{t_1}}(t_1-t_1^{\prime })\cdot e_{l_{t_2}}(t_2-t_2^{\prime }).$$

$$v^i(t)={\left[{\int }_{\mathbb{D}}{g}_1^{i1}(t-s_1)u_1(s_1)d{s}_1\right]}^2+{\left[{\int }_{\mathbb{D}}{g}_1^{i2}(t-s_2)u_1(s_2)d{s}_2\right]}^2,$$

$$h_2^i(t_1;t_2)=g_1^{i1}(t_1)g_1^{i1}(t_2)+g_1^{i2}(t_1)g_1^{i2}(t_2),$$

$$v^i(t)={\int }_{{\mathbb{D}}^2}{h}_2^i(t-s_1;t-s_2)u_1(s_1)u_1(s_2)d{s}_{1}d{s}_2,$$

$${\int }_{t_k^i}^{t_{k+1}^i}{v}^i(t)dt={\kappa }^i{\delta }^i-b^i(t_{k+1}^i-t_k^i).$$

$${\mathcal{T}}_k^i{u}_2=q_k^i,k\in {\mathbb{I}}^i,i=1,\cdots ,M,$$

$${\mathcal{T}}_k^i{u}_2={\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{h}_2^i(t-s_1;t-s_2)u_2(s_1;s_2)d{s}_{1}d{s}_{2}dt,$$

$$\mbox{\boldmath$\Xi$}\mathbf{d}=\mathbf{q},~{}~{}~{}\text{where}~{}~{}~{}\mbox{\boldmath$\Xi$}=[(\mbox{\boldmath$\Xi$}^{1})^{T},...,(\mbox{\boldmath$\Xi$}^{M})^{T}]^{T}~{}~{}~{}\text{and}~{}~{}~{}\mathbf{q}=[(\mathbf{q}^{1})^{T},...,(\mathbf{q}^{M})^{T}]^{T}$$

$$\left[\mbox{\boldmath$\Xi$}^{i}\right]_{k;l_{t_{1}}l_{t_{2}}}=\int_{t^{i}_{k}}^{t^{i}_{k+1}}\!\!\!e_{l_{t_{1}}+l_{t_{2}}}(t)dt\int_{\mathbb{D}^{2}}h^{i}_{2}(s_{1};s_{2})e_{-l_{t_{1}}}(s_{1})e_{-l_{t_{2}}}(s_{2})ds_{1}ds_{2}.$$

$${\hat{u}}_2(t_1;t_2)=\underset{u_2\in {\mathcal{H}}_2}{\mathrm{argmin}}\sum _{i=1}^M\sum _{k\in {\mathbb{I}}^i}({\mathcal{T}}_k^i{u}_2-q_k^i{)}^2.$$

$$\widehat{u_2}(t_1;t_2)=\sum _{l_{t_1}=-L_t}^{L_t}\sum _{l_{t_2}=-L_t}^{L_t}{\hat{d}}_{l_{t_1}{l}_{t_2}}{e}_{l_{t_1}}(t_1)\cdot e_{l_{t_2}}(t_2),$$

$$\hat{\mathbf{d}}=\mbox{\boldmath$\Xi$}^{\dagger}\mathbf{q}$$

$$({\mathcal{P}}_1{h}_1)(t)={\int }_{\mathbb{D}}{h}_1(t^{\prime })K_1(t;t^{\prime })d{t}^{\prime }$$

$$({\mathcal{P}}_2{h}_2)(t_1;t_2)={\int }_{{\mathbb{D}}^2}{h}_2(t_1^{\prime };t_2^{\prime })K_2(t_1,t_2;t_1^{\prime },t_2^{\prime })d{t}_1^{\prime }d{t}_2^{\prime }$$

$${\mathcal{L}}_k^i({\mathcal{P}}_2{h}_2)=q_k^i,k\in {\mathbb{I}}^i,i=1,\cdots ,M,$$

$${\mathcal{L}}_k^i({\mathcal{P}}_2{h}_2)={\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{u}_2^i(t-s_1;t-s_2)(\mathcal{P}{h}_2)(t-s_1;t-s_2)d{s}_{1}d{s}_{2}dt,$$

$$q_k^i={\kappa }^i{\delta }^i-b^i(t_{k+1}^i-t_k^i).$$

$${\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{h}_2(t-s_1;t-s_2)u_2^i(s_1;s_2)d{s}_{1}d{s}_{2}dt=q_k^i.$$

$${\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{h}_2(t-s_1;t-s_2)({\mathcal{P}}_2{u}_2^i)(s_1;s_2)d{s}_{1}d{s}_{2}dt=q_k^i\text{or}$$

$${\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{\int }_{{\mathbb{D}}^2}{h}_2(t-s_1;t-s_2)K_2(s_1,s_2;s_1^{\prime },s_2^{\prime })u_2^i(s_1^{\prime };s_2^{\prime })d{s}_1^{\prime }d{s}_2^{\prime }d{s}_{1}d{s}_{2}dt=q_k^i\text{or}$$

$${\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}{\int }_{{\mathbb{D}}^2}{h}_2(t-s_1;t-s_2)K_2(t-s_1,t-s_2;t-s_1^{\prime },t-s_2^{\prime })d{s}_{1}d{s}_2{u}_2^i(s_1^{\prime };s_2^{\prime })d{s}_1^{\prime }d{s}_2^{\prime }dt=q_k^i\text{or}$$

$${\int }_{t_k^i}^{t_{k+1}^i}{\int }_{{\mathbb{D}}^2}({\mathcal{P}}_2{h}_2)(t-s_1;t-s_2)u_2^i(s_1;s_2)d{s}_{1}d{s}_{2}dt=q_k^i.$$

$${\mathcal{L}}_k^i({\mathcal{P}}_2{h}_2)=q_k^i,k\in {\mathbb{I}}^i,i=1,\cdots ,M.$$

$${\mathcal{P}}_2{h}_2(t_1;t_2)=\sum _{l_{t_1}=-L_t}^{L_t}\sum _{l_{t_2}=-L_t}^{L_t}{h}_{l_{t_1}{l}_{t_2}}{e}_{l_{t_1}}(t_1)\cdot e_{l_{t_2}}(t_2).$$

$$\mbox{\boldmath$\Theta$}\mathbf{h}=\mathbf{q},$$

$$\left[\mbox{\boldmath$\Theta$}^{i}\right]_{k;l_{t_{1}}l_{t_{2}}}=\int_{t^{i}_{k}}^{t^{i}_{k+1}}\!\!\!e_{l_{t_{1}}+l_{t_{2}}}(t)dt\int_{\mathbb{D}^{2}}u^{i}_{2}(s_{1};s_{2})e_{-l_{t_{1}}}(s_{1})e_{-l_{t_{2}}}(s_{2})ds_{1}ds_{2}.$$

$$\widehat{{\mathcal{P}}_2{h}_2}=\underset{{\mathcal{P}}_2{h}_2\in {\mathcal{H}}_2}{\mathrm{argmin}}\sum _{i=1}^M\sum _{k\in {\mathbb{I}}^i}{\left({\mathcal{L}}_k^i({\mathcal{P}}_2{h}_2)-q_k^i\right)}^2.$$

$$\widehat{{\mathcal{P}}_2{h}_2}(t_1;t_2)=\sum _{l_{t_1}=-L_t}^{L_t}\sum _{l_{t_2}=-L_t}^{L_t}{\hat{h}}_{l_{t_1}{l}_{t_2}}{e}_{l_{t_1}}(t_1)\cdot e_{l_{t_2}}(t_2),$$

Full text

Keywords: encoding of visual stimuli, complex cells, quadratic receptive fields, dendritic stimulus processors, sparse neural decoding, sparse functional identification, duality between decoding and functional identification

1. 1 Introduction

2. 2 Neural Circuits with Complex Cells: Encoding, Decoding and Functional Identification

3. 2.1 Modeling Temporal Stimuli

4. 2.2 Encoding of Temporal Stimuli by a Population of Complex Cells

5. 2.3 Decoding of Temporal Stimuli Encoded by a Population of Complex Cells

6. 2.4 Functional Identification of DSPs of Complex Cells

7. 3 Low-Rank Decoding and Functional Identification

8. 3.1 Low-Rank Decoding of Stimuli

9. 3.1.1 Exploiting the Structure of Complex Cell Encoding

10. 3.1.2 Reconstruction Algorithms

11. 3.1.3 Example - Decoding of Temporal Stimuli Encoded with a Population of Complex Cells

12. 3.2 Low-Rank Functional Identification of Complex Cells

13. 3.2.1 Duality Between Low-Rank Functional Identification and Decoding

14. 3.2.2 Functional Identification Algorithms

15. 3.2.3 Example - Identification of Complex Cell DSPs from Spike Times

16. 3.3 Evaluation of Functional Identification of a Neural Circuit of Complex Cells by Decoding

17. 4 Low-Rank Decoding and Functional Identification of Complex Cells with Spatio-Temporal Stimuli

18. 4.1 Modeling of Spatio-Temporal Stimuli

19. 4.2 Encoding of Spatiotemporal Stimuli with a Population of Complex Cells

20. 4.3 Low-Rank Decoding of Spatio-Temporal Visual Stimuli

21. 4.3.1 Example - Decoding of 2D Spatio-Temporal Stimuli

22. 4.3.2 Example - Decoding of 3D Spatio-Temporal Stimuli

23. 4.4 Low-Rank Functional Identification of Spatio-Temporal Complex Cells

24. 4.4.1 Example - Low-Rank Functional Identification of Complex Cell DSP from Spike Times in Response to Spatio-Temporal Stimuli

25. 4.4.2 Example - Evaluation of Functional Identification of Neural Circuit of Complex Cells Using Decoding

26. 4.4.3 Comparison with STC, GQM and NIM

27. 5 Conclusions

1 Introduction

It is widely accepted that the early mammalian visual system employs a series of neural circuits to extract elementary visual features, such as edges and motion [1, 2]. Feature extraction capabilities of simple and complex cells arising in the primary visual cortex (V1) have been extensively investigated. Each simple cell consists of a linear receptive field cascaded with a highly-nonlinear spike generator. Similar to simple cells, complex cells in V1 are selective to oriented edges/lines over a spatially restricted region of the visual field [1]. While simple cells respond maximally to a particular phase of the edge, complex cells are largely phase invariant [3]. Therefore, the receptive fields of complex cells cannot be simply mapped into excitatory and inhibitory regions [1]. Receptive fields of simple cells can be modeled as linear Gabor filters while processing in complex cells can be modeled with a quadrature pair of Gabor filters followed by squaring [4]. Neural circuits comprising complex cells constitute a highly nonlinear circuit as illustrated in Figure 1.

Under the modeling framework of Time Encoding Machines (TEMs) [5, 6, 7], it has been shown that decoding of stimuli and functional identification of linear receptive fields of simple cells are dual to each other [8, 9]. This led to mathematically rigorous identification algorithms for identifying linear receptive fields of simple cells [10]. By modeling the nonlinear processing in complex cells as Volterra Dendritic Stimulus Processors (DSPs) [11, 12], the representation of stimuli encoded by spike times generated by neural circuits with complex cells was also exhaustively analyzed. Functional identification of a complex cell DSP was possible again thanks to the demonstrated duality between decoding and functional identification. While these theoretical methods exhibit deep structural properties, they have been shown to be tractable only for decoding and functional identification problems of small dimensions. In their current form they are not tractable due to the “curse of dimensionality” [13].

The non-linear transformations taking place in the DSP of complex cells lead to loss of phase information. With this in mind, we formulate the reconstruction of stimuli encoded with complex cells as a phase retrieval problem [14] and, in search of tractable algorithms, utilize recent developments in optimization theory of low-rank matrices [15, 16, 14]. By applying such methods, we develop algorithms that are highly effective in decoding visual stimuli encoded by complex cells. As will be detailed in the next sections, the complex cells, as defined in this paper, have DSP kernels that are low-rank and include the ones shown in Figure 1 as a particular case.

After demonstrating that the decoding of visual stimuli becomes tractable, we propose sparse algorithms for functionally identifying the DSPs of complex cells using the spike times they generate. The sparse identification algorithms are based on the key observation that functional identification can be viewed as the dual problem of decoding stimuli that are encoded by an ensemble of complex cells. While the generalization of the duality results from simple cells to complex cells was already given in [11], we show in this paper that these results remain valid under the assumption of sparsity, that is, for the case of low-rank DSP kernels. This significantly reduces the time of stimulus presentation that is needed in the identification process. The sparse duality result also enables us to evaluate the identified circuits in the input space. We achieve the latter by computing the mean square error or signal-to-noise (SNR) of novel stimuli decoded using the identified circuits [8, 9].

This paper is organized as follows. In Section 2, we first introduce the modeling of encoding of temporal stimuli with complex cells. We provide a detailed review of decoding of stimuli and the functional identification of complex cells, and point out the current algorithmic limitations. In Section 3, we provide sparse decoding algorithms that achieve high accuracy and are algorithmically tractable. We then explicate the dual relationship between sparse functional identification and decoding and provide examples for the identification of low-rank, temporal DSP kernels of complex cells. In Section 4, we extend sparse decoding methodology to spatio-temporal stimuli and functional identification of spatio-temporal complex cells. Using novel stimuli, we provide evaluation examples of the identification algorithms in the input space as well as comparisons to other state-of-the-art methods. Finally, we conclude in Section 5 and suggest how the approach advanced in this paper can be applied beyond complex cells.

2 Neural Circuits with Complex Cells:

Encoding, Decoding and Functional Identification

In this section, we model the encoding of temporal stimuli by a neural circuit consisting of neurons akin to complex cells. We start by modeling the space of temporal stimuli in Section 2.1. In Section 2.2, the model of encoding is formally described. In Section 2.3, we proceed to present a reconstruction algorithm for decoding temporal stimuli encoded by the neural circuit. A method for functional identification of neurons constituting the neural circuit is provided in Section 2.4. The reconstruction algorithm and the functional identification algorithm discussed in this section are based on [11].

2.1 Modeling Temporal Stimuli

We model the temporal varying stimuli $u_{1}=u_{1}(t)$, $t\in\mathbb{D}$, to be real-valued elements of the space of trigonometric polynomials [6]. The choice of the space of the trigonometric polynomials has, as we will see, substantial computational advantages.

Definition 1**.**

The space of trigonometric polynomials ${\mathcal{H}}_{1}$ is the Hilbert space of complex-valued functions

[TABLE]

over the domain $\mathbb{D}=[0,S_{t}]$, where

[TABLE]

Here $\Omega_{t}$ denotes the bandwidth, and $L_{t}$ the order of the space. Stimuli $u_{1}\in{\mathcal{H}}_{1}$ are extended to be periodic over ${\mathbb{R}}$ with period $S_{t}=\frac{2\pi L_{t}}{\Omega_{t}}$.

${\mathcal{H}}_{1}$ is a Reproducing Kernel Hilbert Space (RKHS) [17] with reproducing kernel (RK)

[TABLE]

We denote the dimension of ${\mathcal{H}}_{1}$ by $dim({\mathcal{H}}_{1})$ and $dim({\mathcal{H}}_{1})=2L_{t}+1$.

Definition 2**.**

The tensor product space ${\mathcal{H}}_{2}={\mathcal{H}}_{1}\otimes{\mathcal{H}}_{1}$ is a Hilbert space of complex-valued functions

[TABLE]

over the domain $\mathbb{D}^{2}$.

${\mathcal{H}}_{2}$ is an RKHS with reproducing kernel

[TABLE]

Note that $dim({\mathcal{H}}_{2})=dim({\mathcal{H}}_{1})^{2}$.

2.2 Encoding of Temporal Stimuli by a Population of Complex Cells

We consider a neural circuit consisting of $M$ neurons as shown in Figure 2a. For the $i^{\text{th}}$ neuron, input stimulus $u_{1}(t)$ ($u_{1}\in{\mathcal{H}}_{1}$) is first processed by two linear filters with impulse responses $g^{i1}_{1}(t)$ and $g^{i2}_{1}(t)$, the outputs of which are individually squared and then summed together. These processing elements are integral part of the DSP of neuron $i$ [11, 12]. The output of the DSP $i$, denoted by $v^{i}(t)$, is then fed into the Biological Spike Generator (BSG) of neuron $i$. The BSG $i$ encodes the output of DSP $i$ into the spike train $(t^{i}_{k})_{k\in\mathbb{I}^{i}}$. Here $\mathbb{I}^{i}$ is the spike train index set of neuron $i$. We notice the similarity between the overall structure of neural circuits in Figure 2a and Figure 1. In what follows, we refer to the neurons in the neural circuit in Figure 2a as complex cells.

The output of the DSP of the $i^{th}$ neuron in Figure 2a amounts to

[TABLE]

for all $i=1,2,\cdots,M$.

With

[TABLE]

(4) can be rewritten as

[TABLE]

where $\mathbb{D}^{2}$ denotes the Cartesian product of the domain $\mathbb{D}$ of $u_{1}$ and $h^{i}_{2}(t_{1};t_{2})$ is interpreted as a second-order Volterra kernel [18]. We assume that $h_{2}^{i}(t_{1};t_{2})$ is real, bounded-input bounded-output (BIBO) stable, causal and of finite memory. The I/O of the neural circuit shown in Figure 2a can be equivalently outlined as in Figure 2b, in which each neuron processes the input $u_{1}(t)$ nonlinearly by a second order kernel $h^{i}_{2}(t_{1};t_{2})$ followed by a BSG.

Remark 1**.**

Note that the BSG models the spike generation mechanism of the axon hillock of a biological neuron, whereas the DSP is an equivalent model of processing of the stimuli by a sophisticated neural network that proceeds the spike generation. Therefore, stimulus processing and the spike generation mechanism are naturally separated in the neuron model considered here.

For simplicity, we will use the Integrate-and-Fire (IAF) neuron model as the spike generation mechanism (see, e.g., [5]). Note that, the algorithms described here can also be employed with other spike generators such as the Hodgkin-Huxley, Morris-Lecar and Izhikevic neuron models [19, 20, 21, 22, 12]. The integration constant, bias and threshold of the IAF neuron $i=1,2,\cdots,M$, are denoted by $\kappa^{i}$, $b^{i}$ and $\delta^{i}$, respectively. The t-transform of the $i$-th IAF neuron is given by [5, 6, 7]

[TABLE]

Lemma 1**.**

The encoding of the temporal stimulus $u_{1}\in{\mathcal{H}}_{1}$ into the spike train sequence $(t_{k}^{i}),k\in\mathbb{I}^{i}$, $i=1,2,...,M$, by a neural circuit with complex cells is given in functional form by

[TABLE]

where $M$ is the total number of neurons, $n_{i}+1$ is the number of spikes generated by neuron $i$ and $\mathcal{T}^{i}_{k}:{\mathcal{H}}_{2}\rightarrow{\mathbb{R}}$, are bounded linear functionals defined by

[TABLE]

with $u_{2}(t_{1};t_{2})=u_{1}(t_{1})u_{1}(t_{2})$. Finally, $q^{i}_{k}=\kappa^{i}\delta^{i}-b^{i}(t^{i}_{k+1}-t^{i}_{k})$.

Proof: The relationship (8) follows by replacing the functional form of $v^{i}(t)$ given in (6) in equation (7) above. $\Box$

Remark 2**.**

$u_{2}(t_{1},t_{2})=u_{1}(t_{1})\cdot u_{1}(t_{2})$* can be interpreted as a nonlinear map of the stimulus $u_{1}$ into $u_{2}$ defined in a higher dimensional space. The operation performed by the second order Volterra kernel on $u_{2}$ in (9) is linear. Thus, (8) shows that the encoding of temporal stimuli can be viewed as generalized sampling [11].*

2.3 Decoding of Temporal Stimuli Encoded by a Population of Complex Cells

Assuming that the spike times $(t_{k}^{i}),k\in\mathbb{I}^{i},i=1,2,...,M$, are known, by Lemma 1, the neural circuit with complex cells encodes the stimulus via a set of linear functionals acting on $u_{2}$ (see equation (8)). Thus, the reconstruction of $u_{2}$ can in principle be obtained by inverting the set of linear equations (8) [11].

Theorem 1**.**

The coefficients of $u_{2}\in{\mathcal{H}}_{2}$ in (3) satisfy the following system of linear equations

[TABLE]

with $\left[\mathbf{q}^{i}\right]_{k}=q^{i}_{k},\left[\mathbf{d}\right]_{l_{t_{1}}l_{t_{2}}}=d_{l_{t_{1}}l_{t_{2}}}$ and

[TABLE]

The above result can be obtained by plugging (3) into (8). We refer readers to Theorem 1 in [11] for a detailed proof.

We formulate the reconstruction of $u_{2}$ as the following optimization problem:

[TABLE]

Algorithm 1**.**

The solution to (11) is given by

[TABLE]

where $\hat{\mathbf{d}}=[\hat{d}_{{-L_{t}},{-L_{t}}},\cdots,\hat{d}_{{-L_{t}},{L_{t}}},\cdots,\cdots,\hat{d}_{{L_{t}},{-L_{t}}},\cdots,\hat{d}_{{L_{t}},{L_{t}}}]^{T}$ is obtained by

[TABLE]

with † denoting the pseudoinverse operator.

We note that a necessary condition for perfect recovery is that the total number of spikes exceeds $dim({\mathcal{H}}_{1})(dim({\mathcal{H}}_{1})+1)/2+M$ [12]. Therefore, the complexity of the decoding algorithm is on the order of $dim({\mathcal{H}}_{1})^{2}$.

Following [11, 12], the decoding algorithm is called a Volterra Time Decoding Machine (Volterra TDM).

2.4 Functional Identification of DSPs of Complex Cells

In this section, we formulate the functional identification of a single complex cell in the neural circuit described in Figure 2a. We perform $M$ experimental trials. In trial $i,i=1,\cdots,M$, we present a controlled stimulus $u^{i}_{1}(t)$ to the cell and observe the spike times $(t^{i}_{k})_{k\in\mathbb{I}^{i}}$. We assume the cell has a DSP of the from $h_{2}(t_{1};t_{2})=g^{1}_{1}(t_{1})g^{1}_{1}(t_{2})+g^{2}_{1}(t_{1})g^{2}_{1}(t_{2})$ and an integrate and fire BSG with integration constant, bias and threshold denoted by $\kappa,b\mbox{ and }\delta$, respectively. The objective is to functionally identify $h_{2}$ from the knowledge of $u_{1}^{i}$ and the observed spikes $(t^{i}_{k})_{k\in\mathbb{I}^{i}}$, $i=1,\cdots,M$. This is a standard practice in neurophysiology for inferring the functional form of a component of a sensory system [1].

Definition 3**.**

Let $h_{p}\in\mathbb{L}^{1}(\mathbb{D}^{p}),p=1,2$, where $\mathbb{L}^{1}$ denotes the space of Lebesgue integrable functions. The operator $\mathcal{P}_{1}:\mathbb{L}_{1}(\mathbb{D})\rightarrow{\mathcal{H}}_{1}$ given by

[TABLE]

is called the projection operator from $\mathbb{L}^{1}(\mathbb{D})$ to ${\mathcal{H}}_{1}$. Similarly, the operator $\mathcal{P}_{2}:\mathbb{L}_{1}(\mathbb{D}^{2})\rightarrow{\mathcal{H}}_{2}$ given by

[TABLE]

is called the projection operator from $\mathbb{L}^{1}(\mathbb{D}^{2})$ to ${\mathcal{H}}_{2}$.

Note, that for $u_{1}^{i}\in{\mathcal{H}}_{1},\mathcal{P}_{1}u_{1}^{i}=u_{1}^{i}$. Moreover, with $u_{2}^{i}(t_{1},t_{2})=u_{1}^{i}(t_{1})u_{1}^{i}(t_{2}),\mathcal{P}_{2}u_{2}^{i}=u_{2}^{i}$.

Lemma 2**.**

With $M$ trials of stimuli $u^{i}_{2}(t_{1};t_{2})=u^{i}_{1}(t_{1})u^{i}_{1}(t_{2}),i=1,\cdots,M$, presented to a complex cell having DSP $h_{2}(t_{1},t_{2})$, we have

[TABLE]

where

[TABLE]

and

[TABLE]

Proof: With (18) the t-transform for the $i^{th}$ stimulus is given by

[TABLE]

Since $\mathcal{P}_{2}u_{2}^{i}=u_{2}^{i}$, we have

[TABLE]

Finally, with (17), we obtain

[TABLE]

$\Box$

Remark 3**.**

The similarity between equations (8) and (19) suggests that the identification of a complex cell DSP by presenting multiple stimuli is dual to decoding a stimulus encoded by a population of complex cells. This duality is schematically shown in Figure 3.

Theorem 2**.**

Let $\mathcal{P}_{2}h_{2}\in{\mathcal{H}}_{2}$ be of the form

[TABLE]

Then, $\left[\mathbf{h}\right]_{l_{t_{1}}l_{t_{2}}}=h_{l_{t_{1}}l_{t_{2}}}$ with $l_{t_{1}}=-L_{t},\cdots,L_{t},l_{t_{2}}=-L_{t},\cdots,L_{t}$, satisfies the following system of linear equations

[TABLE]

where \mbox{\boldmath\Theta}=[(\mbox{\boldmath\Theta}^{1})^{T},...,(\mbox{\boldmath\Theta}^{M})^{T}]^{T} and $\mathbf{q}=[(\mathbf{q}^{1})^{T},...,(\mathbf{q}^{M})^{T}]^{T}$ with $\left[\mathbf{q}^{i}\right]_{k}=q^{i}_{k},$ and

[TABLE]

Thus, to identify $\mathcal{P}_{2}h_{2}$, we can follow the same methodology as in Algorithm 1, and formulate the functional identification of $\mathcal{P}_{2}h_{2}$ as

[TABLE]

Algorithm 2**.**

The solution to (23) is given by

[TABLE]

where $\hat{\mathbf{h}}=[\hat{h}_{{-L_{t}},{-L_{t}}},\cdots,\hat{h}_{{-L_{t}},{L_{t}}},\cdots,\cdots,\hat{h}_{{L_{t}},{-L_{t}}},\cdots,\hat{h}_{{L_{t}},{L_{t}}}]^{T}$ is obtained by

[TABLE]

The methodology described in Algorithm 25 to identify the nonlinear DSP is called the Volterra Channel Identification Machine (Volterra CIM) [11, 12].

Remark 4**.**

Formulating the decoding and identification problems in the tensor product space ${\mathcal{H}}_{2}$ allows the identification of nonlinear processing by solving a set of linear equations. However, the increased dimensionality necessitates the use of $\mathcal{O}\left(dim\left({\mathcal{H}}_{1}\right)^{2}\right)$ measurements.

3 Low-Rank Decoding and Functional Identification

As shown in Section 2.3, a reconstruction of the signal $u_{2}$ is in principle possible by solving a set of linear equations. However, the complexity of the algorithm is prohibitive. We show in this section that an efficient decoding algorithm can be constructed that exploits the structure of encoding circuits with complex cells. Based on the duality between decoding and functional identification, functional identification algorithms that exploit the structure of the DSP of complex cells are presented. These algorithms largely reduce the complexity of decoding of temporal stimuli encoded by an ensemble of complex cells and that of functional identification of their DSPs.

3.1 Low-Rank Decoding of Stimuli

3.1.1 Exploiting the Structure of Complex Cell Encoding

In Theorem 1, we introduced a vector notation for the coefficients of $u_{2}$

[TABLE]

We introduce here the matrix notation of the coefficients for $u_{2}\in{\mathcal{H}}_{2}$,

[TABLE]

We notice the following: i) since $u_{2}$ is assumed to be real, $\overline{d_{l_{t_{1}},l_{t_{2}}}}=d_{-l_{t_{1}},-l_{t_{2}}}$, and ii) since $u_{2}(t_{1};t_{2})=u_{1}(t_{1})u_{1}(t_{2})=u_{1}(t_{2})u_{1}(t_{1})=u_{2}(t_{2};t_{1})$, we have $d_{l_{t_{1}},l_{t_{2}}}=d_{l_{t_{2}},l_{t_{1}}}$. These properties imply that $\mathbf{D}$ is a Hermitian matrix. Moreover, we note that $u_{2}$ in (8) is the “outer” product of the stimuli $u_{1}$, i.e.,

[TABLE]

where

[TABLE]

are the coefficients of the basis functions of $u_{1}$. Therefore, $\mathbf{D}$ is a rank-1 Hermitian positive semidefinite matrix. This property will be exploited in stimulus decoding (reconstruction).

Theorem 3**.**

Encoding the stimulus $u_{1}\in{\mathcal{H}}_{1}$ with the neural circuit with complex cells given in (7) into the spike train sequence $(t_{k}^{i}),k\in\mathbb{I}^{i}$, $i=1,2,...,M$, satisfies the set of equations

[TABLE]

where $\mbox{\bf Tr}(\cdot)$ is the trace operator, $\mathbf{D}$ is the rank-$1$ positive semidefinite Hermitian matrix $\mathbf{D}=\mathbf{c}\mathbf{c}^{H}$, $q^{i}_{k}=\kappa^{i}\delta^{i}-b^{i}(t^{i}_{k+1}-t^{i}_{k})$ and (\mbox{\boldmath\Phi}^{i}_{k}),k\in\mathbb{I}^{i},i=1,\cdots,M, are Hermitian matrices with entries in the $\left(l_{t_{2}}+L_{t}+1\right)$-th row and $\left(l_{t_{1}}+L_{t}+1\right)$-th column given by

[TABLE]

Proof: Plugging in the general form of $u_{2}$ in (3) into (9), the left hand side of (8) amounts to

[TABLE]

It is easy to verify that the above expression can be written as

[TABLE]

Finally, we note that since $h^{i}_{2},i=1,\cdots,M$, are assumed to be real valued, (\mbox{\boldmath\Phi}^{i}_{k}),k\in\mathbb{I}^{i},i=1,\cdots,M, are Hermitian. $\Box$

Remark 5**.**

We note that equation (30) in Theorem 31 and equation (10) in Theorem 1 are the same. These equations represent the t-transform of a complex cell in (rank-1) matrix and vector form, respectively. The (rank-1) matrix representation is made possible by the equality $u_{2}(t_{1};t_{2})=u_{1}(t_{1})u_{1}(t_{2})$.

3.1.2 Reconstruction Algorithms

Solving the systems of equations (30) and (10) requires at least $dim({\mathcal{H}}_{1})(dim({\mathcal{H}}_{1})+1)/2+M$ measurements. Consequently, practical solutions become quickly intractable. Fortunately, the encoded stimulus is of the form $u_{2}(t_{1};t_{2})=u_{1}(t_{1})u_{2}(t_{2})$. This guarantees that $\mathbf{D}$ is a rank-1 matrix and thus the reconstructed stimulus belongs to a small subset of ${\mathcal{H}}_{2}$. Therefore, we can cast the problem of reconstructing temporal stimuli encoded by neural circuits with complex cells as a feasibility problem, that is, find all positive semidefinite Hermitian matrices that satisfy (30) and have rank 1. As we shall demonstrate, the latter condition can be satisfied with substantially fewer measurements.

Recently, there is an increasing interest in low-rank optimizations such as matrix factorization, matrix completion and rank minimization, both from a theoretical and from a practical standpoint [23, 16, 24]. For example, rank minimization has recently been applied to phase retrieval problems [14].

Our objective here is to find rank-1, positive-semidefinite matrices that satisfy the t-transform (30). Since there always exists at least one rank-1 solution, this is equivalent to the following optimization problem [25]

[TABLE]

The rank minimization problem in (33) is NP-hard. A well known heuristic is to relax the problem (33) to a trace minimization problem [24]. That is, instead of solving (33), we reconstruct $u_{2}$ using Algorithm 36.

Algorithm 3**.**

The reconstruction of $u_{2}$ from the spike times generated by the neural circuit with complex cells is given by

[TABLE]

where

[TABLE]

is the solution to the semidefinite programming (SDP) problem

[TABLE]

When the matrices (\mbox{\boldmath\Phi}^{i}_{k}), $k\in\mathbb{I}^{i},i=1,\cdots,M$, satisfy the rank restricted isometry property [16], the trace norm relaxation converges to the true solution of (33) provided that the number of measurements is of the order \mathcal{O}\Big{(}dim(\mathcal{H}_{1})log\big{(}dim\left(\mathcal{H}_{1}\right)\big{)}\Big{)} [16]. These results suggest that stimuli encoded by complex cells can be decoded with a significantly lower number of measurements than that required by Algorithm 1. To investigate this further, we applied the above algorithm to decode a large number of stimuli encoded by complex cells while varying the number of measurements (spikes) used by the decoding algorithm. The results show that the number of spikes required to faithfully represent a stimulus by a neural circuits consisting of complex cells is quasilinearly rather than quadratically proportional to the dimension of the stimulus space. These results are presented in the subsequent sections.

The matrix of weights $\hat{\mathbf{D}}$ obtained from the above algorithm can be further decomposed to extract the signal $u_{1}$ (up to a sign) as follows.

(i) Perform the eigen-decomposition of $\hat{\mathbf{D}}$. Denote the largest eigenvalue by $\lambda$ and the corresponding eigenvector by $\mathbf{v}$. If (36) does not exactly return a rank-$1$ matrix, choose the largest eigenvalue and disregard the rest. Let $\mathbf{w}=\sqrt{\lambda}\mathbf{v}$.

(ii) The reconstructed stimulus $\hat{u}_{1}$ is given by (up to a sign)

[TABLE]

where

[TABLE]

with $\hat{\mathbf{c}}=\left[\hat{c}_{-L_{t}},\cdots,\hat{c}_{L_{t}}\right]^{T}$, and $[\mathbf{w}]_{L_{t}+1}$ is the $(L_{t}+1)^{\text{th}}$ entry of $\mathbf{w}$, which corresponds to the coefficient $\hat{c}_{0}$.

If $\hat{\mathbf{D}}$ is rank 1, step (i) decomposes $\hat{\mathbf{D}}$ as an “outer” product of a vector and itself (see (28)). The resulting vector $\mathbf{w}$ differs from the actual coefficient vector of the stimulus $u_{1}$ by up to a complex-valued scaling factor. This factor is corrected in step (ii). Since $u_{1}$ is assumed to be real-valued, the “DC” component must be real-valued. Therefore, we rotate $\mathbf{w}$ to remove any imaginary part. In practice, this also ensures $\hat{c}_{-l_{t}}=\overline{\hat{c}_{l_{t}}}$.

Remark 6**.**

Note that we can reconstruct $u_{1}(t)$ up to a sign, since $\mathbf{D}=\mathbf{c}\mathbf{c}^{H}$ and $\mathbf{D}=(-\mathbf{c})(-\mathbf{c}^{H})$ are equally possible.

Remark 7**.**

Note that (33) can be alternatively solved by replacing the objective with the log-det heuristic [24], that is

[TABLE]

where $\lambda>0$ is a small regularization constant. This optimization may further reduce the rank of $\hat{\mathbf{D}}$ when Algorithm 36 fails to progress to an exact rank-1 solution. [24].

While the SDP in (36) provides an elegant way for relaxing the rank minimization problem, it is limited in practice by the need of large amounts of computer memory for numerical calculations. The optimization problem (33) can also be solved using an alternating minimization scheme [26] as outlined in Algorithm 4 below. The alternating minimization approach is more tractable when the dimension of the space is very large. Algorithm 4 uses an initialization step (step $1$ below) that provides an initial iterate whose distance from $\mathbf{D}$ is bounded. It then alternately solves for the left and right singular vector of the rank-1 matrix $\mathbf{D}$ while keeping the other one fixed (step $2$ below). The resulting subproblems admit a straightforward least squares solution, that can be much more efficiently solved than the SDP in Algorithm 36. Moreover, the algorithm is amenable to parallel computation using General Purpose Graphics Processing Units (GPGPUs). The latter property makes it even more attractive when the dimension of the stimulus space is large.

Algorithm 4**.**

Initialize $\hat{\mathbf{c}}_{1}$ and $\hat{\mathbf{c}}_{2}$ to top left and right singular vector respectively of $\sum_{i=1}^{M}\sum_{k\in\mathbb{I}^{i}}q_{k}^{i}\mathbf{\Phi_{k}^{i}}$ normalized to $\sqrt{\frac{1}{\sigma}\sum_{i=1}^{M}\sum_{k\in\mathbb{I}^{i}}(q_{k}^{i})^{2}}$, where $\sigma$ is the top singular value of $\sum_{i=1}^{M}\sum_{k\in\mathbb{I}^{i}}q_{k}^{i}\mathbf{\Phi_{k}^{i}}$. 2. 2.

Solve alternately the following two minimization problems

1. (a)

solve for $\hat{\mathbf{c}}_{1}$ by fixing $\hat{\mathbf{c}}_{2}$

[TABLE] 2. (b)

solve for $\hat{\mathbf{c}}_{2}$ by fixing $\hat{\mathbf{c}}_{1}$

[TABLE]

until \sum_{i=1}^{M}\sum_{\in\mathbb{I}^{i}}(\mbox{\bf Tr}(\mbox{\boldmath\Phi}^{i}_{k}\hat{\mathbf{c}}_{1}\hat{\mathbf{c}}_{2}^{H})-q^{i}_{k})^{2}\leq\epsilon, where $\epsilon>0$ is the error tolerance level. 3. 3.

compute $\hat{\mathbf{D}}=\hat{\mathbf{c}}_{1}\hat{\mathbf{c}}_{2}^{H}$.

$\hat{\mathbf{D}}$ approximates the coefficients of $u_{2}\in{\mathcal{H}}_{2}$ as in (34) We can reconstruct $u_{1}$, using the (appropriately scaled) top eigenvector of $\frac{1}{2}(\hat{\mathbf{D}}+\hat{\mathbf{D}}^{H})$. This can be obtained directly from $\hat{\mathbf{c}}_{1}$ and $\hat{\mathbf{c}}_{2}$ as follows. Let

[TABLE]

and

[TABLE]

the reconstructed stimulus $\hat{u}_{1}$ is given by (up to a sign)

[TABLE]

where

[TABLE]

with $\hat{\mathbf{c}}=\left[\hat{c}_{-L_{t}},\cdots,\hat{c}_{L_{t}}\right]^{T}$.

We point out that we made the decoding manageable by exploiting the structure of $u_{2}$. Therefore, there is no constraint on the exact form $h^{i}_{2}(t_{1};t_{2})$ can take, and the decoding algorithms can be applied to neural circuits with neurons whose DSPs take the form of any second-order Volterra kernel.

3.1.3 Example - Decoding of Temporal Stimuli Encoded with a Population of Complex Cells

Here, the neural circuit we consider consists of $19$ complex cells. The DSPs of the complex cells are of the form

[TABLE]

where $g^{i1}_{1}(t)$ and $g^{i2}_{1}(t)$ are quadrature pairs of temporal Gabor filters and $i=1,\cdots,19$. The Gabor filters are constructed from dilations and translations of the mother wavelets, where the mother functions can expressed as

[TABLE]

and

[TABLE]

The BSG of the complex cells are point IAF neurons with bias $b^{i}=2$ and integration constant $\kappa^{i}=1$, for $i=1,\cdots,M$. These two parameters are kept the same for all stimuli. Different threshold values are chosen for the IAF neurons in order to vary the total number of spikes, which can be used to evaluate how many measurements are required for perfectly reconstructing the input stimuli.

The domain of the input space ${\mathcal{H}}_{1}$ is $\mathbb{D}=[0,1]$ (sec) and $L_{t}=20,\Omega_{t}=20\cdot 2\pi$ (rad/sec). Thus, we have $dim({\mathcal{H}}_{1})=41$. The stimuli were generated by randomly choosing their basis coefficients from an i.i.d. Gaussian distribution.

We tested the encoding and subsequent decoding of $6,570$ stimuli. The total number of spikes produced for each stimulus ranged from 20 to 220. Reconstructions of the stimuli were performed using Algorithm 36, and the SDPs were solved using SDPT3 [27].

We show the SNR of all reconstructions in the scatter plot of Figure 4a. Here solid dots represent exact rank 1 solutions (largest eigenvalue is at least 100 times larger than the sum of the rest of the eigenvalues), and crosses indicate that the trace minimization found a higher rank solution that has a smaller trace. The percentage of exact rank 1 solutions is shown in Figure 4b. A relatively sharp transition from very low probability of recovery to very high rate of perfect reconstruction can be seen, similar to phase transition phenomena in other sparse recovery algorithms [28]. It can also be seen that the number of measurements that are needed for perfect recovery is substantially lower than the $861$ spikes required by decoding based on Theorem 1.

3.2 Low-Rank Functional Identification of Complex Cells

3.2.1 Duality Between Low-Rank Functional Identification and Decoding

As discussed in Section 2.4, the complexity of identification using Algorithm 25 can be prohibitively high. Often, a very large number of stimulus presentation trials are required to fully identify the DSP of biological neurons. To mitigate this, we consider exploiting the structure of the DSP of complex cells as motivated by the tractability of the low rank decoding algorithm when the structure of the stimuli is explored.

We consider a single complex cell whose DSP is of the form

[TABLE]

where $g^{n}_{1}(t),n=1,\cdots,N$, are impulse responses of linear filters, and $N\ll dim({\mathcal{H}}_{1})$. We note that a complex cell described in Figure 2a is a special case of (47) with $N=2$. A natural question here is that, by assuming such a structure, whether the functional identification of complex cell DSPs is tractable.

Remark 8**.**

It is well known that a second-order Volterra kernel has infinite equivalent forms but has a unique symmetric form [18].

We have shown that the low-rank structure of $u_{2}$ leads to a reduction of complexity in the reconstruction of temporal stimuli encoded by an ensemble of complex cells. We also described the duality between decoding and functional identification. If we can show that the functional identification formalism for complex cell DSP is the dual to decoding of low-rank stimuli, it is straightforward to provide tractable algorithms for identifying $h_{2}(t_{1};t_{2})$ of the form (47).

Since $\mathcal{P}_{1}g^{n}_{1}(t)\in{\mathcal{H}}_{1},n=1,\cdots,N$, there is a set of coefficients $(g^{n}_{l_{t}}),l_{t}=-L_{t},...,L_{t}$ and $n=1,2,...,N$, such that

[TABLE]

In what follows wte denote coefficients in vector form as

[TABLE]

Similarly, we denote the coefficients of $\mathcal{P}_{1}h_{2}(t_{1};t_{2})$ in (20) in matrix form as

[TABLE]

Then

[TABLE]

and thus $\mathbf{H}$ is a Hermitian positive semidefinite matrix with rank at most $N$.

Theorem 4**.**

By presenting $M$ trials of stimuli $u^{i}_{2}(t_{1};t_{2})=u^{i}_{1}(t_{1})u^{i}_{1}(t_{2}),i=1,\cdots,M$ to a complex cell, its coefficients satisfy the set of equations

[TABLE]

where $n_{i}+1,i=1,\cdots,M$, is the number of spikes generate by the complex cell in trial $i$, $\mathbf{H}$ is a Hermitian positive semidefinite matrix with $\operatorname{rank}(\mathbf{H})\leq N$, given by $\mathbf{H}=\sum_{n=1}^{N}\mathbf{g}^{n}(\mathbf{g}^{n})^{H},$ with $\mathbf{g}^{n}=\left[g^{n}_{-L_{t}},\cdots,g^{n}_{L_{t}}\right]^{T},$ (\mbox{\boldmath\Psi}^{i}_{k}),k\in\mathbb{I}^{i},i=1,\cdots,M, are Hermitian matrices with entry at $\left(l_{t_{2}}+L_{t}+1\right)$-th row and $\left(l_{t_{1}}+L_{t}+1\right)$-th column given by

[TABLE]

Proof: From Lemma 18, we have

[TABLE]

where

[TABLE]

(52) can be obtained following the steps of the proof of Theorem 31. $\Box$

Remark 9**.**

As in Section 3.2, we note that the similarity in (52) and (30) indicates the duality between low-rank functional identification of complex cells and low-rank decoding of stimuli encoded by a population of complex cells. The duality is illustrated in Figure 5.

3.2.2 Functional Identification Algorithms

To functionally identify the complex cell DSP, we again employ a rank minimization problem

[TABLE]

Algorithm 36 provides a solution to the above rank minimization problem. However, in this case, the optimal solution shall have rank $N$. We relax the problem to a trace minimization problem and consider the following algorithm for low-rank functional identification of complex cells.

Algorithm 5**.**

The functional identification of complex cell DSP from the spike times generated by the neuron in $M$ stimulus trials is given by

[TABLE]

where

[TABLE]

is the solution to the SDP problem

[TABLE]

Based on the results for decoding using Algorithm 36 and provided that $h_{2}$ is of the form (47), we intuitively inferred that the number of measurements for the perfect identification of $\mathcal{P}_{2}h_{2}$ is much smaller than \mathcal{O}\big{(}dim(\mathcal{H}_{1})^{2}\big{)} . We demonstrate that this is the case for a large number of identification examples in the subsequent sections.

This suggests that even if the dimension of the input space becomes large, the functional identification of the DSP of complex cells is still tractable. This result has critical implication for performing neurobiological experiments to functionally identify complex cells. First, it suggests that a much smaller number of stimulus trials is needed for perfect identification. Second, the total number of spikes/measurements that needs to be recorded can be significantly reduced. Both means the duration of experiment can be shortened.

Remark 10**.**

Note that only the projection of the DSP $h_{2}$ onto the space of input stimuli can be identified.

Remark 11**.**

We can use the largest $N$ eigenvalues and their respective eigenvectors of $\hat{\mathbf{H}}$ to obtain the projection of individual linear filter components $\widehat{\mathcal{P}_{1}g^{n}_{1}},n=1,\cdots,N$. However, these components may not directly correspond to $\mathcal{P}_{1}g^{n}_{1},n=1,\cdots,N$, in that the original projections may not be “orthogonal”, whereas the eigenvalue decomposition imposes orthogonality.

As in Algorithm 4 when applied for solving the decoding problem, the rank minimization problem above can be solved using alternating minimization, as described in Algorithm 6 below. Here, we solve for the top $N$ left and right singular vectors of $\mathbf{H}$ alternately, where $N$ is the rank of the second order Volterra DSP. We note that the initialization step is akin to running an algorithm very similar to the spike-triggered covariance (STC) algorithm widely used in neuroscience [29, 30, 31, 32, 33]. The subsequent steps then improve upon this initial estimate.

Algorithm 6**.**

Initialize $\hat{\mathbf{H}}_{1}$ and $\hat{\mathbf{H}}_{2}$ to top $N$ left and right singular vectors, respectively, of $\sum_{i=1}^{M}\sum_{k=1}^{n_{i}}q_{k}^{i}\mathbf{\Psi_{k}^{i}}$ with the $n^{th}$ singular vector normalized to $\frac{1}{N}\sqrt{\frac{1}{\sigma_{n}}\sum_{i=1}^{M}\sum_{k=1}^{n_{i}}(q_{k}^{i})^{2}}$, where $\sigma_{n}$ is the top $n^{th}$ singular value of $\sum_{i=1}^{M}\sum_{k=1}^{n_{i}}q_{k}^{i}\mathbf{\Psi_{k}^{i}}$. 2. 2.

Solve the following two minimization problems

1. (a)

solve for $\hat{\mathbf{H}}_{1}$ by fixing $\hat{\mathbf{H}}_{2}$

[TABLE] 2. (b)

solve for $\hat{\mathbf{H}}_{2}$ by fixing $\hat{\mathbf{H}}_{1}$

[TABLE]

until \sum_{i=1}^{M}\sum_{\in\mathbb{I}^{i}}(\mbox{\bf Tr}(\mbox{\boldmath\Psi}^{i}_{k}\hat{\mathbf{H}}_{1}\hat{\mathbf{H}}_{2}^{H})-q^{i}_{k})^{2}\leq\epsilon, where $\epsilon>0$ is the error tolerance level. 3. 3.

compute $\hat{\mathbf{H}}=\frac{1}{2}\left(\hat{\mathbf{H}}_{1}\hat{\mathbf{H}}_{2}^{H}+\hat{\mathbf{H}}_{2}\hat{\mathbf{H}}_{1}^{H}\right)$.

3.2.3 Example - Identification of Complex Cell DSPs from Spike Times

In this example, we consider identifying a single complex cell having the following Volterra DSP

[TABLE]

where

[TABLE]

In repeated trials we presented to the complex cell 1-second long stimuli chosen from the input space. The domain of the input space ${\mathcal{H}}^{1}_{1}$ is $\mathbb{D}=[0,1]$ (sec) and $L_{t}=20,\Omega_{t}=20\cdot 2\pi$ (rad/sec) and thus, $dim({\mathcal{H}}^{1}_{1})=41$. The stimuli were generated by independently choosing their basis coefficients from the same Gaussian distribution. We presented a total of $16,600$ different stimuli in the repeated trials. We then randomly selected between 30-80 trial subsets such that the total number of spikes in each subset was between $60$ and $160$. We performed the identification process on each subset using Algorithm 59. The optimization problem was solved using SDPT3.

For each instantiation of the identification algorithm, we recorded whether the optimization process resulted in a rank-2 solution and also the SNR of the identified DSP with respect to the original one. For the purpose of demonstration, we binned these results based on number of spikes used into bins of width $10$. The percentage of rank-2 solutions is shown in Figure 6a as a function of number of measurements. The mean SNR is shown in Figure 6b.

It can be seen from Figure 6b that the identification algorithm presented here is able to recover the underlying DSP with exceptional accuracy using a reasonable and tractable number of measurements.

3.3 Evaluation of Functional Identification of a Neural Circuit of Complex Cells by Decoding

In Section 3.1, we have shown that the sparse decoding algorithm requires much less number of neurons and measurements (spikes) in the reconstruction of stimuli encoded by a neural circuit of complex cells. We have also demonstrated in Section 3.2 that the proposed sparse functional identification algorithm enables the identification of complex cells with a tractable number of measurements. Together, the two algorithms afford us tractable functional identification of an entire neural circuit of complex cells that is capable of fully representing stimuli information, in that i) the size of the neural circuit is tractable, and ii) the requirement for functional identification is tractable.

In [9] and [10], it was shown that the evaluation of functional identification of an entire neural circuit can be more intuitively performed in the input space by decoding the stimuli with identified circuit parameters. Here, we extend the previous results and apply such evaluation procedure on the sparse decoding and sparse functional identification algorithms. The procedure is described as follows. First, each complex cell is functionally identified using Algorithm 59 or Algorithm 6. Second, novel stimuli are presented to the neural circuit. Third, the spike trains observed are used to reconstruct the encoded novel stimuli by the sparse decoding algorithm, assuming that the circuit parameters take the identified values. Finally, SNR of the reconstruction can be obtained. A high SNR indicates a well identified circuit while a low number implies that the functional identification of the neural circuit is not of good quality. The latter can be caused by a lack of number of measurements used in functional identification, or by a lack of complex cells in the neural circuit.

We performed the functional identification of all $19$ complex cells in the neural circuit given in the example in Section 3.1.3. We first identified all complex cells by presenting to the neural circuit $M$ temporal stimuli. We repeated the identification of the entire circuit using $8$ different values of $M$. We then presented to the same circuit (with the original DSPs as in Section 3.1.3), 100 novel stimuli drawn from the input space and used the spike times generated by the neural circuit to decode the stimuli. In the decoding process however, we assumed that the DSPs of the set of complex cells are as identified, for all $8$ values of $M$. The mean reconstruction SNR of the $100$ stimuli is shown in Figure 7. As shown, the quality of reconstruction is low until enough trials were used in identification. When more than 19 trials were performed, perfect reconstruction of the entire neural circuit was achieved. The dimension of the stimulus space was $41$ and the average number of spikes per neuron used for identification varied from $44$ for $6$ trials to $202$ for $28$ trials.

4 Low-Rank Decoding and Functional Identification of Complex Cells with Spatio-Temporal Stimuli

In this section, we extend our results obtained in Section 3 to neural circuits with complex cells that encode spatio-temporal stimuli. We will first introduce the space of spatio-temporal stimuli in Section 4.1. In Section 4.2, we formulate the encoding of spatio-temporal stimuli by a population of complex cells. By extending Theorem 31 to Theorem 5, we argue in Section 4.3 that decoding of stimuli encoded by a neural circuit of complex cells is tractable for spatio-temporal stimuli. We then employ extension of the Algorithms 36 and 4 developed in Section 3.1.2 and demonstrate their effectiveness with a few examples. In Section 4.4, we investigate the duality between functional identification of spatio-temporal DSPs of complex cells and decoding of stimuli encoded by complex cells with a bank of spatio-temporal DSPs. Here we extend Theorem 53 to the encoding circuits with spatio-temporal complex cells, Theorem 6. We then apply extensions of the Algorithms 59 and 6 developed in Section 3.2.2 to the identification of spatio-temporal DSPs of complex cells and demonstrate their effectiveness.

4.1 Modeling of Spatio-Temporal Stimuli

The stimuli $u_{1}$ defined here have $p$ spatial dimensions and a single temporal dimension, i.e., $u_{1}=u_{1}(x_{1},x_{2},\cdots,x_{p},t)$. For simplicity of notation, we use a compact, vector notation and denote the spatial variables as $\mathbf{x}=(x_{1},x_{2},\cdots,x_{p})$. When $p=2$, $u_{1}$ is the usual visual stimulus.

Definition 4**.**

The space of trigonometric polynomials ${\mathcal{H}}_{1}^{p}$ is the Hilbert space of complex-valued functions

[TABLE]

where

[TABLE]

[TABLE]

over the domain $\mathbb{D}$, where, by abuse of notation, $\mathbb{D}=[0,S_{x_{1}}]\times[0,S_{x_{2}}]\times\cdots[0,S_{x_{p}}]\times[0,S_{t}]$ and $S_{t}=\frac{2\pi L_{t}}{\Omega_{t}}\;,S_{x_{1}}=\frac{2\pi L_{x_{1}}}{\Omega_{x_{1}}}\;,S_{x_{2}}=\frac{2\pi L_{x_{2}}}{\Omega_{x_{2}}},\cdots,S_{x_{p}}=\frac{2\pi L_{x_{p}}}{\Omega_{x_{p}}}$. In addition, $e_{\mathbf{l_{x}}l_{t}}(\mathbf{x},t)=e_{\mathbf{l_{x}}}(\mathbf{x})e_{l_{t}}(t)$ where

[TABLE]

and

[TABLE]

Here $\Omega_{t}$ denotes the bandwidth, and $L_{t}$ the order of the space in the temporal domain while $\Omega_{x_{i}}$ and $L_{x_{i}}$ denote the bandwidth and order of the space in the $i^{th}$ spatial variable. Stimuli $u_{1}\in{\mathcal{H}}_{1}^{p}$ are periodic with periods $S_{t},S_{x_{1}},\cdots,S_{x_{p}}$.

${\mathcal{H}}_{1}^{p}$ is a Reproducing Kernel Hilbert Space (RKHS) with reproducing kernel (RK)

[TABLE]

We denote the temporal dimension of ${\mathcal{H}}_{1}^{p}$ by $dim_{t}({\mathcal{H}}_{1}^{p})=2L_{t}+1$ and the total dimension by $dim({\mathcal{H}}_{1}^{p})=(2L_{t}+1)\prod_{i=1}^{p}(2L_{x_{i}}+1)$.

Definition 5**.**

The tensor product space ${\mathcal{H}}_{2}^{p}={\mathcal{H}}_{1}^{p}\otimes{\mathcal{H}}_{1}^{p}$ is a Hilbert space of complex-valued functions

[TABLE]

over the domain $\mathbb{D}^{2}$.

${\mathcal{H}}_{2}^{p}$ is an RKHS with reproducing kernel

[TABLE]

Note that $dim({\mathcal{H}}_{2}^{p})=(dim({\mathcal{H}}_{1}^{p}))^{2}$.

4.2 Encoding of Spatiotemporal Stimuli with a Population of Complex Cells

We consider again a neural circuit consisting of a population of $M$ neurons modeling a population of complex cells as illustrated in Figure 1. The input to the neural circuit is spatiotemporal stimulus as defined in Section 4.1.

The input stimulus $u_{1}(\mathbf{x},t)$ to neuron $i$ is first processed by two spatio-temporal linear filters whose impulse responses are denoted, by abuse of notation, as $g^{i1}_{1}(\mathbf{x},t)$ and $g^{i2}_{1}(\mathbf{x},t)$, respectively. The output of the linear filters are squared and summed. The sum $v^{i}(t)$, as the output of the DSP, is then fed into the BSG of neuron $i$. The BSG encodes the DSP output into the spike train $(t^{i}_{k})_{k\in\mathbb{I}^{i}}$. Here $\mathbb{I}^{i}$ is the spike train index set of neuron $i$.

Similar to the temporal case, the neural circuit is equivalent to that shown in Figure 8a. Here, the output of the DSP for each neuron $i=1,2,\cdots,M$, can be expressed as

[TABLE]

Here

[TABLE]

is the low-rank DSP [11]. The encoding of stimulus by the neural circuit with complex cells is a special case of the low-rank DSP of the form given in (70). When using IAF point neurons as models of the BSGs, we have the following theorem describing the encoding of stimuli.

Lemma 3**.**

The encoding of stimulus $u_{1}\in{\mathcal{H}}^{p}_{1}$ into the spike train sequence $(t_{k}^{i}),k\in\mathbb{I}^{i},i=1,2,...,M,$ by a neural circuit of spatio-temporal complex cells is given in functional form by

[TABLE]

where $\mathcal{T}^{i}_{k}:{\mathcal{H}}_{2}^{p}\rightarrow{\mathbb{R}}$, are bounded linear functionals defined by

[TABLE]

with $u_{2}(\mathbf{x_{1}},t_{1};\mathbf{x_{2}},t_{2})=u_{1}(\mathbf{x_{1}},t_{1})u_{1}(\mathbf{x_{2}},t_{2})$. Finally, $q^{i}_{k}=\kappa^{i}\delta^{i}-b^{i}(t^{i}_{k+1}-t^{i}_{k})$.

Proof: As in Lemma 1, the t-transform of the $i$-th IAF neuron is given by (7).

The relationship (71) follows after replacing $v^{i}(t)$ given in (69) in equation (7). $\Box$

Similar to Remark 2, equation (71) shows that the encoding of a stimuli by the neural circuit with low-rank DSPs can be viewed as generalized sampling.

By abuse of notation, we denote by $\mathbf{c}$ the vector representing the coefficients of $u_{1}$ in (65), and $\mathbf{D}$ as the matrix representing the coefficients of $u_{2}$ in (67). We skip here the detailed entries of $\mathbf{c}$ and $\mathbf{D}$ due to the complexity of the indices, but their construction follows closely with (29) and (27), respectively, and $\mathbf{D}=\mathbf{c}\mathbf{c}^{H}$.

Theorem 5**.**

Encoding the stimulus $u_{1}\in{\mathcal{H}}^{p}_{1}$ with the neural circuit with complex cells given in (69) into the spike train sequence $(t_{k}^{i}),k\in\mathbb{I}^{i}$, $i=1,2,...,M$, satisfies the set of equations

[TABLE]

*where $\mathbf{D}=\mathbf{c}\mathbf{c}^{H}$ is a rank-$1$ Hermitian matrix and (\mbox{\boldmath\Phi}^{i}_{k}), $k\in\mathbb{I}^{i},i=1,\cdots,M$, are Hermitian matrices. [\mbox{\boldmath\Phi}^{i}_{k}]_{\mathbf{l_{x_{2}}}l_{t_{2}};\mathbf{l_{x_{1}}}l_{t_{1}}} denotes the entry at the

$\left((l_{t_{2}}+L_{t_{2}}+1)\prod_{i=1}^{p}(L_{x_{i2}}+1)+\sum_{j=1}^{p}(l_{x_{j2}}+L_{x_{j2}}+1)\prod_{i=1}^{j-1}(2L_{x_{i2}}+1)\right)$-th row and the

$\left((l_{t_{1}}+L_{t_{1}}+1)\prod_{i=1}^{p}(L_{x_{i1}}+1)+\sum_{j=1}^{p}(l_{x_{j1}}+L_{x_{j1}}+1)\prod_{i=1}^{j-1}(2L_{x_{i1}}+1)\right)$-th column, and*

[TABLE]

where $\mathbf{l_{x_{i}}}=(l_{x_{1i}},l_{x_{2i}},\cdots,l_{x_{pi}}),i=1,2$.

Proof: Plugging in the general form of $u_{2}$ in (67) into (72), the left hand side of (71) amounts to

[TABLE]

It is easy to verify that the expression above can be written as

[TABLE]

where the

$\left((l_{t_{1}}+L_{t_{1}}+1)\prod_{i=1}^{p}(L_{x_{i1}}+1)+\sum_{j=1}^{p}(l_{x_{j1}}+L_{x_{j1}}+1)\prod_{i=1}^{j-1}(2L_{x_{i1}}+1)\right)$-th row

$\left((l_{t_{2}}+L_{t_{2}}+1)\prod_{i=1}^{p}(L_{x_{i2}}+1)+\sum_{j=1}^{p}(l_{x_{j2}}+L_{x_{j2}}+1)\prod_{i=1}^{j-1}(2L_{x_{i2}}+1)\right)$-th column entry of $\mathbf{D}$ amounts to $\left[\mathbf{D}\right]_{\mathbf{l_{x_{1}}}l_{t_{1}};\mathbf{l_{x_{2}}}l_{t_{2}}}=d_{\mathbf{l_{x_{1}}},l_{t_{1}},-\mathbf{l_{x_{2}}},-l_{t_{2}}}$.

Since $u_{2}(\mathbf{x_{1}},t_{1};\mathbf{x_{2}},t_{2})=u_{1}(\mathbf{x_{1}},t_{1})u_{1}(\mathbf{x_{2}},t_{2})$ and $d_{\mathbf{l_{x_{1}}},l_{t_{1}},-\mathbf{l_{x_{2}}},-l_{t_{2}}}=c_{\mathbf{l_{x_{1}}},l_{t_{1}}}c_{\mathbf{l_{x_{2}}},l_{t_{2}}}^{H}$, thereby $\mathbf{D}=\mathbf{c}\mathbf{c}^{H}$. We also note that since $h^{i}_{2},i=1,\cdots,M$, are assumed to be real valued, (\mbox{\boldmath\Phi}^{i}_{k}),k\in\mathbb{I}^{i},i=1,\cdots,M, are Hermitian. $\Box$

4.3 Low-Rank Decoding of Spatio-Temporal Visual Stimuli

When using an algorithm similar to Algorithm 1 to reconstruct spatio-temporal stimuli encoded by a neural circuit with complex cells, at least $dim({\mathcal{H}}^{p}_{1})\left(dim({\mathcal{H}}^{p}_{1})+1\right)/2$ measurements are required. In addition, at least $dim({\mathcal{H}}^{p}_{1})\left(dim({\mathcal{H}}^{p}_{1})+1\right)/(4L_{t}+1)$ neurons are required, a number that can become unrealistically high with an increasing dimension of the input space.

With the observation that $\mathbf{D}=\mathbf{cc}^{H}$ is a rank-one matrix, we can apply algorithms similar to those described in Section 3.1.2 to recover spatio-temporal stimuli encoded by a population of spiking neurons with low-rank DSPs. For the sake of brevity, we skip the details of the extended Algorithms 36 and 4, and in what follows we will provide some examples that demonstrate that the decoding of spatio-temporal stimuli is still tractable.

4.3.1 Example - Decoding of 2D Spatio-Temporal Stimuli

We first present an example in which $\mathbf{x}$ is one-dimensional, i.e., $\mathbf{x}=x_{1}$. In this example, our main focus is to illustrate how the number of spikes affects the reconstruction of stimuli encoded by complex cells.

The neural circuit we consider here consists of 62 direction selective complex cells. The low-rank DSPs of the complex cells are of the form

[TABLE]

where $g^{i1}_{1}(x_{1},t)$ and $g^{i2}_{1}(x_{1},t)$ are quadrature pairs of spatio-temporal Gabor filters and $i=1,\cdots,M$. The Gabor filters are constructed from dilations and translations of the mother wavelets on a dyadic grid, where the mother functions can expressed as

[TABLE]

and

[TABLE]

The BSG of the complex cells are IAF neurons with bias $b^{i}=10$ and integration constant $\kappa=1$, for $i=1,\cdots,M$. These two parameters are kept the same for all stimuli. Different threshold values are chosen for the IAF neurons in order to vary the total number of spikes in a larger range to evaluate how many measurements are required for a perfect reconstruction of input stimuli.

The domain of the input space ${\mathcal{H}}^{1}_{1}$ is $\mathbb{D}=[0,32]\times[0,0.4]$ ([a.u.] and [sec], respectively) and $L_{x_{1}}=6,L_{t}=4,\Omega_{x_{1}}=0.1875\cdot 2\pi,\Omega_{t}=10\cdot 2\pi$ [rad/sec]. Thus, $dim({\mathcal{H}}^{1}_{1})=117$. Stimuli were randomly generated by choosing the basis coefficients to be i.i.d. Gaussian random variables.

We tested the encoding of $1,416$ stimuli. Each time, a different number of spikes was generated. The reconstruction of stimuli was performed in MATLAB using the extended Algorithm 36, and the SDPs were solved using SDPT3 [27].

The SNR of all reconstructions is depicted in the scatter plot of Figure 9a. Here solid dots represent exact rank 1 solutions (largest eigenvalue is at least 100 times larger than the sum of the rest of the eigenvalues), and crosses indicate that the trace minimization found a higher rank solution with a smaller trace. The percentage of exact rank 1 solutions is shown in Figure 9b. Similar to phase transition phenomena in other sparse recovery algorithms [28], a relatively sharp transition (around 50 spikes) from very low probability of recovery to very high probability of perfect reconstruction can be seen. It can also be seen that the number of measurements that are needed for perfect recovery is substantially lower than the $6,965$ spikes required by Algorithm 1.

4.3.2 Example - Decoding of 3D Spatio-Temporal Stimuli

Next, we present two examples of decoding of spatio-temporal visual stimuli encoded by a population of complex cells. Here, $\mathbf{x}=(x_{1},x_{2})$ and the Volterra DSPs of the complex cells are of the form

[TABLE]

where $g^{i1}_{1}(\mathbf{x},t)$ and $g^{i2}_{1}(\mathbf{x},t)$ are, for simplicity, quadrature pairs of spatial-only Gabor filters and $i=1,\cdots,M$. The Gabor filters are constructed from dilations, translations and rotations of a mother wavelets [6],

[TABLE]

and

[TABLE]

For the first example, a 0.4-second-long synthetically generated video sequence is encoded by the neural circuit. The order of the input space was chosen to be $L_{x_{1}}=L_{x_{2}}=3,L_{t}=4$. Thus, the dimension of the input space is $441$. The input stimulus was created by choosing its basis coefficients to be i.i.d. Gaussian random variables. The stimulus was encoded by a neural circuit consists of $318$ complex cells. A total of $1,374$ spikes were generated by the encoding circuit. The stimulus was decoded using the extended Algorithm 36. As shown in Figure 10, the video sequence can be perfectly reconstructed with a fairly small number of spikes (A snapshot of the video is shown, see also Supplementary Video S1 for full video). The SNR of the reconstructed video was $92.8$ [dB], thereby reaching almost perfect reconstruction with machine precision. Note that without the reconstruction algorithm employed here, $97,461$ measurements would be required from at least $5,733$ complex cells to achieve perfect reconstruction.

The second example uses a natural video sequence. As an illustration, we project the natural video into a space with lower spatial bandwidth in order to reduce the dimension of the embedding space. Here, the order of the space is given by $L_{x_{1}}=6,L_{x_{2}}=6,L_{t}=2$, and thus, the dimension of input space is $845$. A neural circuit consisting of $472$ complex cells was used to encode the stimulus, and a total of $6,000$ spikes were generated. The stimulus was reconstructed using the extended Algorithm 4. The number of spikes employed, $\sim 6,000$, is much lower than the $397,150$ spikes required by Algorithm 1 for perfect recovery. A snapshot of the original and reconstructed video sequence are shown in Figure 11 (see also Supplementary Video S2). The SNR of the reconstructed video was $68.0$ [dB].

Both examples demonstrate the effectiveness of the reconstruction algorithm proposed in this paper. In particular, compared to the decoding algorithm based on Theorem 1, the number of neurons in the neural circuit was substantially reduced. Therefore, with biologically plausible spike rates and number of complex cells, the information contained in the spike times of these neurons can faithfully represent visual stimuli.

4.4 Low-Rank Functional Identification of Spatio-Temporal Complex Cells

Similar to Section 3.2, we consider here the identification of low-rank DSP of complex cells from spike times generated when multiple stimulus trials are presented. We first define the projection operators in ${\mathcal{H}}^{p}_{1}$. Then, based on (69), we show that the duality between decoding and functional identification also holds in the spatio-temporal case.

Definition 6**.**

Let $h_{n}\in\mathbb{L}^{1}(\mathbb{D}^{n}),n=1,2$, where $\mathbb{L}^{1}$ denotes the space of Lebesgue integrable functions. The operator $\mathcal{P}_{1}^{p}:\mathbb{L}_{1}(\mathbb{D})\rightarrow{\mathcal{H}}^{p}_{1}$ given by

[TABLE]

is called the projection operator from $\mathbb{L}^{1}(\mathbb{D})$ to ${\mathcal{H}}^{p}_{1}$. Similarly, the operator $\mathcal{P}_{2}^{p}:\mathbb{L}_{1}(\mathbb{D}^{2})\rightarrow{\mathcal{H}}_{2}$ given by

[TABLE]

is called the projection operator from $\mathbb{L}^{1}(\mathbb{D}^{2})$ to ${\mathcal{H}}_{2}$.

We consider here complex cells whose low-rank DSP can be expressed more generally as

[TABLE]

where, by abuse of notation, $g^{n}_{1}(\mathbf{x},t),n=1,\cdots,N$ are impulse responses of spatio-temporal linear filters, and $N\ll dim({\mathcal{H}}^{p}_{1})$. Similar to the approach we take in Section 3.2, this particular structure can be exploited to identify the projection of $h_{2}$ using tractable algorithms.

By abuse of notation, we denote $\mathbf{g}^{n}$ as the vector representing the coefficients of $\mathcal{P}^{p}_{1}g^{n}_{1}$, and $\mathbf{H}$ as the matrix representing the coefficients of $\mathcal{P}^{p}_{2}h_{2}$. The detailed entries of $\mathbf{g}^{n}$ and $\mathbf{H}$ are constructed similarly to (49) and (50), respectively. In addition, we have $\mathbf{H}=\sum_{n=1}^{N}\mathbf{g}^{n}(\mathbf{g}^{n})^{H}$.

Theorem 6**.**

By presenting $M$ trials with stimuli $u^{i}_{2}(\mathbf{x_{1}},t_{1};\mathbf{x_{2}},t_{2})=u^{i}_{1}(\mathbf{x_{1}},t_{1})u^{i}_{1}(\mathbf{x_{2}},t_{2}),i=1,\cdots,M$, to a complex cell and observing the spike trains $t^{i}_{k},k\in\mathbb{I}^{i},i=1,2,\cdots,M$, the coefficients of the projections $\mathcal{P}_{2}^{p}h_{2}$ of the DSP of the complex cell, satisfy the set of equations

[TABLE]

*where $\mathbf{H}$ is a rank-$N$ positive semidefinite Hermitian matrix and (\mbox{\boldmath\Psi}^{i}_{k}), $k\in\mathbb{I}^{i},i=1,\cdots,M$, are Hermitian matrices with the entry at the

$\left((l_{t_{2}}+L_{t_{2}}+1)\prod_{i=1}^{p}(L_{x_{i2}}+1)+\sum_{j=1}^{p}(l_{x_{j2}}+L_{x_{j2}}+1)\prod_{i=1}^{j-1}(2L_{x_{i2}}+1)\right)$-th row and the

$\left((l_{t_{1}}+L_{t_{1}}+1)\prod_{i=1}^{p}(L_{x_{i1}}+1)+\sum_{j=1}^{p}(l_{x_{j1}}+L_{x_{j1}}+1)\prod_{i=1}^{j-1}(2L_{x_{i1}}+1)\right)$-th column given by [\mbox{\boldmath\Psi}^{i}_{k}]_{\mathbf{l_{x_{2}}}l_{t_{2}};\mathbf{l_{x_{1}}}l_{t_{1}}}=*

[TABLE]

where $\mathbf{l_{x_{i}}}=(l_{x_{1i}},l_{x_{2i}},\cdots,l_{x_{pi}}),i=1,2$.

Proof: Essentially similar to the proof of Theorem 53.

Remark 12**.**

Theorem 5 and Theorem 6 suggest that decoding of spatio-temporal stimuli encoded by a population of complex cells is dual to the functional identification of the DSP of complex cells presented with multiple stimulus trials. This is further illustrated in Figure 8. Note that in identification, only the projection of the complex cell DSP onto the stimulus space can be identified.

Based on Theorem 6, we can provide functional identification algorithms for complex cell DSPs of the form (84) with a significant reduction in the number of required trials and spikes. The algorithms are similar to those presented in Section 3.2.2. In what follows we present a few example of identification of DSPs of complex cells.

4.4.1 Example - Low-Rank Functional Identification of Complex Cell DSP from Spike Times in Response to Spatio-Temporal Stimuli

In this example, we first consider identifying the DSP of a single complex cell in the neural circuit used in Section 4.3.1. As a reminder, the neural circuit used in the example in Section 4.3.1 encodes spatio-temporal stimuli of the form $u_{1}(x_{1},t)$.

We presented to the population of $M$ complex cells 0.4-second stimuli, where $M$ varied from $40$ to $80$. The stimuli were generated by choosing their basis coefficients as i.i.d. Gaussian random variables. For each $M$, we repeated the functional identification process for $200$ times, each with different stimuli. Identification was essentially based on the extended Algorithm 36, where the SDPs were again solved by SDPT3.

The percentage of rank 2 solutions is shown in Figure 12a as a function of number of experimental trials. The mean SNR is shown in Figure 12b. Figure 12a suggests that, if the number of trials is larger than $70$, the solution to the trace minimization coincides with high probability with the rank minimization problem. In contrast, identification of the complex cell DSP using Algorithm 25 would have required at least $407$ trials.

It can be easily seen that the identification process does not require a large number of trials to achieve perfect identification, thereby enabling the identification of non-linear dendritic processing of cells similar in structure to complex cells with a tractable amount of physiological recordings.

4.4.2 Example - Evaluation of Functional Identification of Neural Circuit of Complex Cells Using Decoding

We then performed the functional identification of all $62$ complex cells in the neural circuit used of the example in Section 4.3.1. Here, our goal is to evaluate the identification quality using decoding.

We first identified all complex cells by presenting to the neural circuit $M$ spatio-temporal stimuli. We also performed the identification of the entire circuit using $8$ different values of $M$. We then presented to the same circuit 100 novel stimuli drawn from the input space and used the spike times generated by the neural circuit to decode the stimuli. In the decoding process, we assumed that the DSPs of the set of complex cells are as identified, for all $8$ values of $M$. The mean reconstruction SNR of the $100$ stimuli is shown in Figure 13. As shown, the quality of reconstruction was kept at low SNR until enough trials were used in identification. When more than 70 trials were performed, perfect reconstruction was achieved, and thereby the entire neural circuit has been identified with a very high quality.

4.4.3 Comparison with STC, GQM and NIM

We compared the performance of the low-rank functional identification algorithm introduced here with the widely used Spike-Triggered Covariance (STC) algorithm [31]. As in Section 4.4.1, a complex cell with a pair of orthogonal Gabor filters was chosen for identification. However, the filters had different norms.

Figure 14a shows the quality of identification (SNR) as the number of spikes used in identification increases. Note that the low-rank functional identification algorithm reached perfect identification using only 746 spikes, whereas the performance of the STC algorithm saturated at $\sim 17$ [dB] after almost $40,000$ spikes were used. Figure 14b shows the identified individual Gabor filters of the complex cells using both algorithms. The number of spikes used are indicated at the top of each column.

We also evaluated the identification performance of the generalized quadratic model (GQM) [34] and the non-linear input model (NIM) [35] with quadratic upstream filters to the same example above. The results (not shown) were similar to those obtained with the STC algorithm.

We note that while the low-rank functional identification algorithm is formulated as non-linear sampling using TEMs and solved using recent advances in low-rank matrix sensing, the other algorithms tested here rely on moment based or likelihood based methods that require a large number of samples to converge.

5 Conclusions

In this paper, we presented sparse algorithms for the reconstruction of temporal as well as spatio-temporal stimuli from spike times generated by neural circuits consisting of complex cells. We developed these algorithms by exploiting the structure of complex cells with low-rank DSP kernels and shown that the reconstruction algorithms become tractable. For neural circuits consisting of complex cells, this suggests that, in addition to each extracting visual features, a biologically plausible number of complex cells are capable of faithfully representing visual stimuli.

Based on duality between sparse decoding and functional identification, we showed that functional identification of complex cells DSPs can be efficiently achieved using similar algorithms as used in decoding. These algorithms makes the functional identification of complex cells tractable, allowing guaranteed high quality identification using a much smaller set of testing stimuli as well as of shorter time duration.

The mathematical treatment presented here, however, is not limited to the complex cells in V1. It can be applied to other neural circuits of interest. For example, early olfactory coding in fruit flies [36] and auditory encoding in grasshoppers [37] have also been shown to have the structure of low-rank DSP kernels. Moreover, the Hassenstein-Reichardt detector [38], a popular model for elementary motion detectors in fruit flies, is also I/O equivalent to low-rank DSP kernels.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

The research reported here was supported by AFOSR under grant #FA9550-16-1-0410.