Automating the design of informative sequences of sensory stimuli

Abstract

Adaptive stimulus design methods can potentially improve the efficiency of sensory neurophysiology experiments significantly; however, designing optimal stimulus sequences in real time remains a serious technical challenge. Here we describe two approximate methods for generating informative stimulus sequences: the first approach provides a fast method for scoring the informativeness of a batch of specific potential stimulus sequences, while the second method attempts to compute an optimal stimulus distribution from which the experimenter may easily sample. We apply these methods to single-neuron spike train data recorded from the auditory midbrain of zebra finches, and demonstrate that the resulting stimulus sequences do in fact provide more information about neuronal tuning in a shorter amount of time than do more standard experimental designs.

Appendices

Appendix A: Using a frequency representation to smooth the STRF

To represent the STRF in the Fourier domain, we applied the Fourier transform separately to the spectral and temporal dimensions of the STRF. Applying the separable Fourier transform to the STRF is just a linear transformation. This transformation maps the STRF into a coordinate system in which the basis functions are rank one matrices. Each of these matrices is the product of 1-dimensional sine-waves in the spectral and temporal directions of the STRF. Using these basis functions we can write the STRF such that each row and column of the STRF is a linear combination of 1-d sine-waves,

$$\begin{array}{lll} \theta(i,j)&=&\sum\limits_{\alpha=1}^{{\ensuremath{m_f}}}\sum\limits_{\beta=1}^{{\ensuremath{m_t}}} {\ensuremath{\gamma}}^1_{\alpha,\beta}\sin(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i) \\ &&\times\sin(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j) \end{array} $$

(30)

$$\begin{array}{lll} &&+\sum\limits_{\alpha=1}^{{\ensuremath{m_f}}}\sum\limits_{\beta=0}^{{\ensuremath{m_t}}} {\ensuremath{\gamma}}^2_{\alpha,\beta}\sin(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i)\\ &&\times\cos(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j)\\ \end{array} $$

(31)

$$\begin{array}{lll} &&+ \sum\limits_{\alpha=0}^{{\ensuremath{m_f}}}\sum\limits_{\beta=1}^{{\ensuremath{m_t}}} {\ensuremath{\gamma}}^3_{\alpha,\beta}\cos(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i) \\ &&\times\sin(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j) \\ \end{array} $$

(32)

$$\begin{array}{lll} &&+\sum\limits_{\alpha=0}^{{\ensuremath{m_f}}}\sum\limits_{\beta=0}^{{\ensuremath{m_t}}} {\ensuremath{\gamma}}^4_{\alpha,\beta}\cos(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i)\\ &&\times \cos(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j). \end{array} $$

(33)

The functions $\sin(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i)$ and $\cos(2\pi\cdot{\ensuremath{f_{o,f}}}\cdot\alpha\cdot{}i)$ determine how each basis function varies across the spectral dimension of the STRF while the functions $\sin(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j)$ and $\cos(2\pi\cdot{\ensuremath{f_{o,t}}}\cdot\beta\cdot{}j)$ determine how the basis functions vary across time in the STRF. Each pair of sine-waves measures the amount of energy at particular frequencies in the spectral and temporal dimensions. The amplitude of each frequency is determined by the coefficients ${\ensuremath{\gamma}}^i_{{\ensuremath{\alpha}},{\ensuremath{\beta}}}$. To form an orthogonal basis for the STRF we need to project the STRF onto sinusoids with frequencies

$$\{0,{\ensuremath{f_{o,f}}},2{\ensuremath{f_{o,f}}},\ldots {\ensuremath{m_f}}{\ensuremath{f_{o,f}}}\} \quad \{0,{\ensuremath{f_{o,t}}},2{\ensuremath{f_{o,t}}},\ldots,{\ensuremath{m_t}}{\ensuremath{f_{o,t}}}\} $$

(34)

$${\ensuremath{f_{o,f}}}=\frac{1}{{\ensuremath{n_f}}} \quad {\ensuremath{f_{o,t}}}=\frac{1}{{\ensuremath{n_t}}} \\ $$

(35)

$${\ensuremath{m_f}}=\left\lceil\frac{1}{2{\ensuremath{f_{o,f}}}}-1\right\rceil \quad {\ensuremath{m_t}}=\left\lceil\frac{1}{2{\ensuremath{f_{o,t}}}}-1\right\rceil; $$

(36)

${\ensuremath{f_{o,f}}}$ and ${\ensuremath{f_{o,t}}}$ are the fundamental frequencies and are set so that one period corresponds to the dimensions of the STRF (${\ensuremath{n_t}}$ and ${\ensuremath{n_f}}$ denote the dimensions of the STRF in the time and frequency dimensions, respectively), and ${\ensuremath{m_f}}$ and ${\ensuremath{m_t}}$ are the largest integers such that ${\ensuremath{m_f}}{\ensuremath{f_{o,f}}}$ and ${\ensuremath{m_t}}{\ensuremath{f_{o,t}}}$ are less than the Nyquist frequency. We subtract 1 and take the ceiling to make sure the frequencies of our basis functions are less than the Nyquist frequency. The unknown parameters in this new coordinate system are the amplitudes, ${\ensuremath{\vec{\gamma}}}=\{{\ensuremath{\gamma}}^1_{\alpha,\beta},{\ensuremath{\gamma}}^2_{\alpha,\beta}, {\ensuremath{\gamma}}^3_{\alpha,\beta},{\ensuremath{\gamma}}^4_{\alpha,\beta}\}$. For simplicity, we will continue to refer to the unknown parameters as θ, realizing that the STRF is represented using this new basis. Since this transformation is linear we can continue to apply our methods for fitting the GLM and optimizing the stimuli.

To low pass filter the STRF we can simply force the coefficients of ${\ensuremath{\theta}}$ corresponding to high frequencies to zero; i.e we pick cutoffs ${\ensuremath{n_{tc}}}$ and ${\ensuremath{n_{fc}}}$ for the time and spectral directions respectively and set

$${\ensuremath{\gamma}}^i_{{\ensuremath{\alpha}},{\ensuremath{\beta}}}=0 \quad \textrm{i}f \,\, \alpha > {\ensuremath{n_{fc}}} \,\, \textrm{or} \,\, \beta > {\ensuremath{n_{tc}}}. $$

(37)

Decreasing the cutoff frequencies not only makes the estimated STRFs smoother, it also reduces the dimensionality of the model. Reducing the dimensionality makes it easier to fit the GLM and optimize the stimuli, but the risk is that the lower-dimensional model may be too simple to adequately model auditory neurons. We can mitigate this risk by using a soft cutoff. Rather than force all high-frequency components to zero, we can adjust our prior to reflect our strong belief that high frequencies should have little energy; we simply set the prior mean of these coefficients to zero and decrease their prior variance. If we now estimate the STRF using the maximum of the posterior then the amplitudes at high frequencies will be biased by our prior towards zero. However, given sufficient evidence the posterior mean will yield non-zero estimates for the amplitudes of high frequencies. See Theunissen et al. (2001) for details and David et al. (2007), Calabrese et al. (2010) for further discussion.

We chose to impose a hard cutoff because we wanted to reduce the dimensionality to make online estimation of the model and online optimization of the stimuli more tractable. To pick the cutoff frequencies, we picked a single neuron and estimated the STRF using maximum-likelihood for a variety of cutoff frequencies. We evaluated the quality of each model by computing the log-likelihood of the bird’s responses to inputs in a test set. The test set consisted of one bird song and one ml-noise stimulus which were not used to train the models. We chose the cutoff frequencies to be ${\ensuremath{n_{fc}}}=10$ and ${\ensuremath{n_{tc}}}=4$ because these values provided good predictive performance for both the bird song and ml-noise while keeping the number of unknown parameters tractable (in this case the STRF has 189 unknown parameters).

Appendix B: Computing the average information for a Gaussian process

In this section we show how the average information per stimulus,

$${\ensuremath{E_{{\ensuremath{\theta}}}}} \log \left| E_{s} \exp(s^T{\ensuremath{\theta}})ss ^T \right|, $$

can be computed when the input distribution is a Gaussian process. For the GLM the expected Fisher information matrix ${\ensuremath{E_{s}}}J_{\rm exp}(s,{\ensuremath{\theta}})$ has a simple 1-dimensional dependence on ${\ensuremath{\theta}}$,

$$J_{\rm exp}(s,{\ensuremath{\theta}})=J_{\rm exp}(s^T{\ensuremath{\theta}})ss^T $$

(38)

$$J_{\rm exp}(s^T{\ensuremath{\theta}})=-E_{r} \frac{\partial^2 \log p(r|\rho=s^T{\ensuremath{\theta}})}{\partial \rho^2} ss^T \\ $$

(39)

$$ =J_{\rm exp}(\rho=s^T{\ensuremath{\theta}})ss^T . $$

(40)

This 1-dimensional structure along with the fact that p(s) is Gaussian makes computing the expectations tractable. We start by defining a new coordinate system in which the first axis is aligned with ${\ensuremath{\theta}}$. This coordinate system is defined by the orthonormal matrix, ${\mathcal{R}_{\theta}}$. The first column of ${\mathcal{R}_{\theta}}$ is $\frac{{\ensuremath{\theta}}}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}}$ and the remaining columns are a suitable set of orthonormal vectors. We can thus define the transformation of s and ${\ensuremath{\theta}}$ into this new coordinate system,

$${\ensuremath{\theta^{'}}}={\mathcal{R}_{\theta}}^T{\ensuremath{\theta}} $$

(41)

$${\ensuremath{\vec{w}}}={\mathcal{R}_{\theta}}^T s. $$

(42)

This coordinate system has the convenient properties

$${\ensuremath{\theta^{'}_{i}}}=0 \quad \forall i\neq1 \\ $$

(43)

$${\Rightarrow}{\ensuremath{\vec{w}}}{}^T{\ensuremath{\theta^{'}}} = w_{1}{\ensuremath{\theta^{'}_{1}}}.$$

(44)

We can now rewrite our objective function

$${\ensuremath{\mathcal{F}}}(p(s))={\ensuremath{E_{{\ensuremath{\theta}}}}}\log| {\ensuremath{E_{p(s)}}} J_{\rm exp}(\rho)ss^T| $$

(45)

$$= {\ensuremath{E_{{\ensuremath{\theta}}}}}\log| {\ensuremath{E_{p({\ensuremath{\vec{w}}})}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}){\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}{}^T| $$

(46)

$$\label{eqn:optfun} ={\ensuremath{E_{{\ensuremath{\theta}}}}}\log| {\ensuremath{E_{{w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}){\ensuremath{E_{{{w_2}},\ldots,{w_{\rm{dim}(s)}|w_{1}}}}}{\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}{}^T|. $$

(47)

Since p(s) is Gaussian and ${\ensuremath{\vec{w}}}={\mathcal{R}_{\theta}}^Ts$, $p({\ensuremath{\vec{w}}})$ is Gaussian with mean ${\mathcal{R}_{\theta}}^T {\mu_s}$ and covariance matrix ${\mathcal{R}_{\theta}}^T {C_s} {\mathcal{R}_{\theta}}^T$. Consequently, $p({\ensuremath{\vec{w}}}|w_{1})$ is also Gaussian and can be computed using the standard Gaussian conditioning formulas,

$$\begin{array}{lll} p({\ensuremath{\vec{w}}}|w_{1})&=&\mathcal{N}({\mathcal{R}_{\theta}}^T{\mu_s}+\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\mathcal{R}_{\theta}}^T{\ensuremath{\boldsymbol{\gamma}}}(w_{1}-{\ensuremath{\mu_{\omega_1}}}),\\ &&{\mathcal{R}_{\theta}}^T{C_s}{\mathcal{R}_{\theta}}-\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\mathcal{R}_{\theta}}^T{\ensuremath{\boldsymbol{\gamma}}}{\ensuremath{\boldsymbol{\gamma}}}^T{\mathcal{R}_{\theta}}) \end{array} $$

(48)

$${\ensuremath{\mu_{\omega_1}}}=\frac{{\ensuremath{\theta}}^T}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}}{\mu_s} $$

(49)

$${\ensuremath{\sigma^2_{\omega_1}}}=\frac{{\ensuremath{\theta}}^T}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}}{C_s} \frac{{\ensuremath{\theta}}}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}} \\ $$

(50)

$${\ensuremath{\boldsymbol{\gamma}}}= {C_s} \frac{{\ensuremath{\theta}}}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}}. $$

(51)

Using this distribution we can easily compute the conditional expectation,

$$\begin{array}{lll} &&{\kern-12pt}{\ensuremath{E_{{{\ensuremath{\vec{w}}}|w_{1}}}}} {\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}^T \\ &&= {\mathcal{R}_{\theta}}^T\left({C_s} - \frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\boldsymbol{\gamma}}}{\ensuremath{\boldsymbol{\gamma}}}^T+\left({\mu_s}+\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\boldsymbol{\gamma}}}(w_{1}-{\ensuremath{\mu_{\omega_1}}})\right)\vphantom{\left(\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}\right)^T}\right. \\&&{\kern3pc}\quad\times\left.\left({\mu_s}+\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\boldsymbol{\gamma}}}(w_{1}-{\ensuremath{\mu_{\omega_1}}})\right)^T \right){\mathcal{R}_{\theta}}\\ \end{array} $$

(52)

$$\begin{array}{lll} &&={\mathcal{R}_{\theta}}^T\left({C_s} +\left({\mu_s}+\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\boldsymbol{\gamma}}}(w_{1}-{\ensuremath{\mu_{\omega_1}}})-\frac{1}{\sqrt{{\ensuremath{\sigma^2_{\omega_1}}}}}{\ensuremath{\boldsymbol{\gamma}}}\right)\vphantom{\left(\frac{1}{\sqrt{{\ensuremath{\sigma^2_{\omega_1}}}}}\right)^T}\right. \\ &&\quad{\kern3pc}\times\left.\left({\mu_s}+\frac{1}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\boldsymbol{\gamma}}}(w_{1}-{\ensuremath{\mu_{\omega_1}}})+\frac{1}{\sqrt{{\ensuremath{\sigma^2_{\omega_1}}}}}{\ensuremath{\boldsymbol{\gamma}}}\right)^T \right){\mathcal{R}_{\theta}}\\ \\ \end{array} $$

(53)

$$ ={\mathcal{R}_{\theta}}^T\left({C_s} + \left({\ensuremath{\vec{\kappa}}{w}}_{1}+{\ensuremath{\vec{\delta}}}\right)\left({\ensuremath{\vec{\kappa}}{w}}_{1}+{\ensuremath{\vec{\eta}}}\right)^T\right){\mathcal{R}_{\theta}} $$

(54)

$${\ensuremath{\vec{\kappa}}}=\frac{{\ensuremath{\boldsymbol{\gamma}}}}{{\ensuremath{\sigma^2_{\omega_1}}}}\\ $$

(55)

$${\ensuremath{\vec{\delta}}}={\mu_s}-\frac{{\ensuremath{\boldsymbol{\gamma}}}}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\mu_{\omega_1}}}-\frac{{\ensuremath{\boldsymbol{\gamma}}}}{\sqrt{{\ensuremath{\sigma^2_{\omega_1}}}}}\\ $$

(56)

$${\ensuremath{\vec{\eta}}}={\mu_s}-\frac{{\ensuremath{\boldsymbol{\gamma}}}}{{\ensuremath{\sigma^2_{\omega_1}}}}{\ensuremath{\mu_{\omega_1}}}+\frac{{\ensuremath{\boldsymbol{\gamma}}}}{\sqrt{{\ensuremath{\sigma^2_{\omega_1}}}}}\\ $$

(57)

The key point is that the expected value is just a rank-1 perturbation of a rotated C _s . We can now evaluate the expectation over w ₁,

$$\begin{array}{lll} &&{\kern-20pt}{\ensuremath{E_{{w_{1}}}}}J_{\rm exp}\left(w_{1}{\ensuremath{\theta^{'}_{1}}}\right){\ensuremath{E_{{{\ensuremath{\vec{w}}}|w_{1}}}}} {\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}^T\\ &&\;={\mathcal{R}_{\theta}}^T\left( {C_s}{\ensuremath{\varpi_1}}+ {\ensuremath{\varpi_3}}\left[\left({\ensuremath{\vec{\kappa}}}+\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}{\ensuremath{\vec{\delta}}}\right)\left({\ensuremath{\vec{\kappa}}}+\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}{\ensuremath{\vec{\eta}}}\right)^T\right.\right.\\ &&\qquad\qquad\qquad\qquad\qquad\left.\left. +\left(\frac{{\ensuremath{\varpi_1}}}{{\ensuremath{\varpi_3}}}-\left(\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}\right)^2\right){\ensuremath{\vec{\delta}}}{\ensuremath{\vec{\eta}}}^T\right] \right){\mathcal{R}_{\theta}}\\ \label{eqn:ejww} \\ \end{array} $$

(58)

$${\ensuremath{\varpi_1}}={\ensuremath{E_{{w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}})\\ $$

(59)

$${\ensuremath{\varpi_2}}={\ensuremath{E_{{w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}) w_{1} \\ $$

(60)

$${\ensuremath{\varpi_2}}={\ensuremath{E_{{w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}) w_{1}^2; $$

(61)

$w_{1}=\frac{{\ensuremath{\theta}}^T}{{\ensuremath{||{{\ensuremath{\theta}}}||_2}}}s$, and p(w ₁) is Gaussian with mean and variance $({\ensuremath{\mu_{\omega_1}}},{\ensuremath{\sigma^2_{\omega_1}}})$. The above are just 1-dimensional expectations so for any value of ${\ensuremath{\theta}}$ we could compute them numerically.

Equation (59) is a rank 2 update of C _s . Therefore we can use the matrix determinant lemma to compute $|{\ensuremath{E_{{w_{1}}}}}{\ensuremath{E_{{{\ensuremath{\vec{w}}}|w_{1}}}}} J_{\rm exp} {\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}^T|$,

$$\begin{array}{lll} &&{\kern-16pt}\log|{\ensuremath{E_{{w_{1}}}}}{\ensuremath{E_{{{\ensuremath{\vec{w}}}|w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}) {\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}^T| \\ &&=\rm{dim}({C_s})\log {\ensuremath{\varpi_1}} + \log | I + V^T({\ensuremath{\varpi_1}}{C_s})^{-1}U| \\ &&{\kern.8pc} +\log |{C_s}| \end{array} $$

(62)

$$U={\ensuremath{\varpi_3}} \left[\left({\ensuremath{\vec{\kappa}}}+\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}{\ensuremath{\vec{\delta}}}\right), \left(\frac{{\ensuremath{\varpi_1}}}{{\ensuremath{\varpi_3}}}-\left(\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}\right)^2\right){\ensuremath{\vec{\delta}}}\right]\\ $$

(63)

$$V=\left[\left({\ensuremath{\vec{\kappa}}}+\frac{{\ensuremath{\varpi_2}}}{{\ensuremath{\varpi_3}}}{\ensuremath{\vec{\eta}}}\right), {\ensuremath{\vec{\eta}}}\right]. $$

(64)

Since $I + V^T({\ensuremath{\varpi_1}}{C_s})^{-1}U$ is a 2-d matrix, we can compute its determinant analytically. Taking the expectation with respect to ${\ensuremath{\theta}}$ yields,

$$\begin{array}{lll} &&{\kern-6pt}{\ensuremath{E_{{\ensuremath{\theta}}}}}\log|{\ensuremath{E_{{w_{1}}}}}{\ensuremath{E_{{{\ensuremath{\vec{w}}}|w_{1}}}}} J_{\rm exp}(w_{1}{\ensuremath{\theta^{'}_{1}}}) {\ensuremath{\vec{w}}}{\ensuremath{\vec{w}}}^T| \\ &&{\kern6pt}=\rm{dim}({C_s}){\ensuremath{E_{{\ensuremath{\theta}}}}}\log {\ensuremath{\varpi_1}} + {\ensuremath{E_{{\ensuremath{\theta}}}}}\log | I \\ &&{\kern15pt}+ V^T({\ensuremath{\varpi_1}}{C_s})^{-1}U|+\log |{C_s}|. \end{array} $$

(65)

Appendix C: Solving the optimization problem described in Section 4.1.3

Our goal is to solve the optimization problem

$$ {\rm{arg}}\max\limits_{{\mu_s},{C_s}} d{\ensuremath{\vec{\mu}_{t}}}^T{\mu_s}+\frac{d}{2}Tr({C_s}{}R)+ \log|{C_s}|, $$

(66)

where we have abbreviated

$$ R={\ensuremath{\vec{\mu}_{t}}}{\ensuremath{\vec{\mu}_{t}}}^T+C_t, $$

(67)

over all (μ _s ,C _s ) subject to the constraints

$$s^T{C_s}s > 0 \quad \forall s\neq0 $$

(68)

$$Tr({C_s})< m- ||{\mu_s}||^2, $$

(69)

where m is the maximum allowed average stimulus power. Clearly the optimal μ _s will be parallel to ${\ensuremath{\vec{\mu}_{t}}}$. Therefore, only the magnitude of the optimal μ _s is unknown. We can therefore rewrite the objective function as

$${\rm{arg}}\max\limits_{||{\mu_s}||} \left[ d ||{\ensuremath{\vec{\mu}_{t}}}||||{\mu_s}|| +{\rm{arg}}\max_{{C_s}} \left(\frac{d}{2}Tr({C_s}{}R)+ \log|{C_s}| \right) \right]$$

(70)

We rewrite the inner problem to make the dependence on ||μ _s || (through the power constraint) more explicit:

$$ {\rm{arg}}\max\limits_{{C_s}}\frac{d}{2}Tr({C_s}{}R)+ \log|{C_s}| $$

(71)

$$s.t \quad s^T{C_s}s > 0 \quad \forall s\neq0 \\ $$

(72)

$$Tr({C_s})< m- ||{\mu_s}||^2. $$

(73)

We solve this optimization problem by introducing a Lagrange multiplier,

$$ L=\frac{d}{2}Tr({C_s}{}R)+\log|{C_s}| - \lambda Tr({C_s})$$

(74)

$$=Tr \left[ {C_s}{} \left( \frac{d}{2} (R-\lambda{}I ) \right) \right] + \log|{C_s}|.$$

(75)

This Lagrangian is exactly isomorphic to twice the log-likelihood in a multivariate Gaussian model with zero mean, if we interpret C _s as the inverse covariance matrix in the Gaussian model and $- \frac{d}{2} (R-\lambda{}I )$ as the observed sample covariance matrix. Standard arguments involving a change of basis now imply that the optimal C _s is given by

$$ {C_s} = -\left( \frac{d}{2} (R-\lambda{}I ) \right)^{-1}, \label{eq:cs} $$

(76)

for any λ >  max _i {r _i }, where r _i denotes the i-th eigenvalue of R. This condition on λ is required to ensure that the resulting C _s maximizes the Lagrangian L, and guarantees that C _s is positive definite.

We now solve for λ by plugging this C _s into our power constraint:

$$ Tr({C_s})=\frac{2}{d}\sum_i \frac{1}{\lambda-r_i} = m - ||{\mu_s}||^2 \label{eq:lambda} $$

(77)

We can easily solve this equation numerically on the allowed range λ >  max _i r _i to compute λ as a function of ||μ _s ||. We can then in principle do a search over all ||μ _s || to find the optimal value (μ _s ,C _s ). In fact, a more efficient method is to instead just compute the optimal (||μ _s ||,C _s ) for each value of λ; thus, a single 1-d search over λ is guaranteed to find the optimal (μ _s ,C _s ). Also note that we can compute the inverse in Eq. (77) efficiently for any value of λ by computing the eigendecomposition of R once and then using the fact that eig(R − λI) = eig(R) − λ; we used this formula already in Eq. (77).

Enforcing stationarity by incorporating Toeplitz constraints

It is worth noting, in the case of stimulus filters θ that extend over more than one time bin (i.e., n _t as defined in Appendix A is greater than one), that the stimulus sequence drawn from the Gaussian distribution defined above will not be temporally stationary. Instead, the stimulus sequence will consist of a series of appended n _t -long segments of draws from a Gaussian distribution, and therefore the marginal distribution of the inputs s _t will in general be an n _t -mixture of Gaussians, instead of a single Gaussian distribution. We recover a single Gaussian only in the special stationary case that the stimulus covariance C _s is constrained to have a Toeplitz structure and the mean vector μ _s is constrained to be constant with respect to time-shifts.

Since we are observing the neural responses r _t given the input s _t presented at each time point t, we should arguably optimize our information over this marginal mixture distribution, instead of the single Gaussian distribution optimized in this appendix. Alternately, we could enforce stationarity in our optimized Gaussian process by including Toeplitz constraints on C _s in our derivation above. We have had limited success deriving a computationally efficient optimization strategy in either of these cases, but this remains an attractive direction for future research. Meanwhile, the results described in Section 4.2 (with a non-stationary optimized Gaussian stimulus ensemble) remain encouraging.