Choosing parameters when using natural movies

[jsiegle]

I am trying to fit STRFs for visual cortex spike train data collected under conditions of natural movie stimulation.

The stimulus is an N x 45 x 60 movie with a 30 Hz frame rate. I set the response to an array of size N, with a 1 whenever a spike occurred, and 0 otherwise.

There's clearly structure in the spike-triggered average, but so far I haven't been able to get anything useful from splineLG. Do you have a sense for which fitting method and parameters space would work well for this application?

[huangziwei]

Since your response is spike, it would be better to use splineLNP. I haven't got the time to write an example for LNP yet, but it is almost the same as splineLG, but with one extra hyperparameter dt. Here's a quick example:

from rfest import splineLNP
lnp = splineLNP(X, y, dt=0.033, dims=[45, 60], df=15)

When I get new data and have no ideas if there will be good RFs, I will normally play around with df first, and plot w_sta and w_spl to see if there are visible RFs or not. If yes, then add regularizations using lnp.fit().

I never got the chance to play with natural stimuli, so I can't say for sure that it will works for your case. I will keep this issue open.

[jsiegle]

I chose a cell with a well-defined response in the STA, but I'm having trouble getting anything out of the spline fit, regardless of the df setting. Here's a typical example:

Any suggestions for other things to try?

[huangziwei]

Interesting. But I am not sure if your STA is well defined.

I'd be careful when interpreting STA estimated from natural movies, since the power spectrum of your stimulus is not flat (meaning pixels in movies are highly correlated), and the STA will likely be distorted. But of course, the quality of STA depends on many other factors, such as the shape of the RF (if it is Gaussian-like, it will likely be bad; but if it's Gabar-like, will probably be good). See my simulation here: https://gist.github.com/huangziwei/de026bbc2fbe3535a14cf3444e909df3

From my experience, sometimes choosing a wrong number of time lag will also corrupt the estimated RF. Since your stimulus is quite fast (30 Hz), maybe try increasing it to -0.5s or -0.8s or even -1s?

[philippberens]

@jsiegle is that the raw STA without prewhitening?

[huangziwei]

I also encountered cases where w_spl is crap while w_opt looks fine (and this happened more often with splineLNP). So you can also try runninglnp.fit(num_iters=200, alpha=1, lambd=0.05, verbal=20) and also play with lambd to see if it will work.

This simulation (assuming softplus nonlinearity and poisson noise) looks similar to your results:

[jsiegle]

Thanks for the advice, I think I'm getting closer.

First of all, for all cells the STA is usually dominated by the average movie frame. The input frames are prewhitened, but most of the STAs end up looking like this (10 frames prior to the spike are shown):

It's only when you take the derivative of the STA that you can see structure emerging:

I used the cell above as input to the splineLNP object, and after running lnp.fit, it seems like it's returning something reasonable:

This is using df=8 and lambd=0.05. Do you have any suggestions for further improving things, or does this look good to you?

[huangziwei]

Cool. I think all you need to do now is to try different combinations of those hyperparameters. Your RF looks Gabor-ish, it normally requires larger df than Gaussian-like RFs to fit its shape.

In splineLNP, you can also set nonlinearity='exponential' (default is softplus). Some papers say exponential function works better than softplus, but for me it's generally the other way around. But it never hurt to try.

[jsiegle]

I'm still looking into hyperparameter tuning, but I ran into a problem when I switched to using GPU computation: the design matrix won't fit in GPU RAM unless I use many fewer frames, or spatially downsample. Have you encountered this issue? Is there any workaround besides switching back to the CPU-only version of jax? I have 128 GB of CPU RAM, but only 8 GB on the GPU.

[huangziwei]

I mainly use the CPU-only version, and only do some tests to see if it runs on GPU or not. So I don't have much experience on that.

But I am wondering at which step the memory error occurred. Previously I set compute_mle=True by default. This is the most expensive step, which eats up quite a lot of memory and not so useful most of the time. Maybe you can try setting it to False. Or just pull the latest push, I just changed the default.

[jsiegle]

Ok, good to know, I will use the CPU-only version instead.

It crashes when building the splineLNP object if X is larger than 8 GB, presumably when it tries to create a DeviceArray of this size on the GPU. For 10 time lags, this happens if the stimulus is longer than ~40k frames (out of 162k in the full dataset).