fail to estimate dpss_windows for long signals about nitime HOT 14 CLOSED

Hi Alex,

From looking at the error:

In [12]: alg.dpss_windows(166800,4,8)

ValueError Traceback (most recent call last)

/Users/arokem/ in ()

/Library/Frameworks/Python.framework/Versions/6.0.0/lib/python2.6/site-packages/nitime/algorithms/spectral.py
in dpss_windows(N, NW, Kmax)
420 ab[1] = ((N - 1 - 2 * nidx) / 2.) ** 2 * np.cos(2 * np.pi * W)
421 # only calculate the highest Kmax-1 eigenvectors

--> 422 l, v = linalg.eig_banded(ab, select='i', select_range=(N -
Kmax, N - 1))
423 dpss = v.transpose()[::-1]
424

/Library/Frameworks/Python.framework/Versions/6.0.0/lib/python2.6/site-packages/scipy/linalg/decomp.pyc
in eig_banded(a_band, lower, eigvals_only, overwrite_a_band, select,
select_range, max_ev)
529 range = select, lower = lower,
530 overwrite_ab = overwrite_a_band,
--> 531 abstol=abstol)
532 # crop off w and v

533         w = w[:m]

ValueError: array is too big.

I am guessing that the lapack function that is doing the heavy lifting
in there (that's 'bevx' on line 531 of the scipy.linalg file) is
trying to make an N by N matrix and that is too large for numpy(?):

In [16]: n = 166800

In [17]: np.empty((n,n))

ValueError Traceback (most recent call last)

/Users/arokem/ in ()

ValueError: array is too big.

I am not sure how to solve this, though.

On Thu, May 5, 2011 at 8:20 AM, agramfort
[email protected]
wrote:

I have times series with 166800 samples (raw MEG data).

alg.dpss_windows(166800, 4, 8)

fails. However, in Matlab it works.

any idea of how to fix this?

Reply to this email directly or view it on GitHub:
#61

Ariel Rokem
Helen Wills Neuroscience Institute
University of California, Berkeley
http://argentum.ucbso.berkeley.edu/ariel

from nitime.

yes the pb occurs in the lapack function. The pb is therefore relatively deep, hence will require
a better understanding the numerical tricks that can be used to compute the tapers. In Matlab
the key functions seems to be tridieig. It might require to go deep in the scipy lapack bindings...

from nitime.

Alright, I am back from the conference and I have been looking around a bit for clues on this. It turns out that we are not the only ones worried about the efficiency of finding eigenvalues of large tri-diagonal matrices with eig_banded:

http://mail.scipy.org/pipermail/scipy-user/2011-February/028550.html

With a little bit of creative googling, I managed to find code taken from this book:

http://www.amazon.com/gp/product/0521191327/ref=s9_simh_gw_p14_d0_i1?pf_rd_m=ATVPDKIKX0DER&pf_rd_s=center-2&pf_rd_r=1SJ9T50258Z23E78DRCE&pf_rd_t=101&pf_rd_p=470938631&pf_rd_i=507846

Which seems to do something like that. I put the code here:

https://gist.github.com/970864

I am not sure what license this under and also, it's pure python, which makes me suspect that it might not be as good as wrapping the relevant lapack functions mentioned in that thread. I haven't tried it out yet, so I don't even know whether it does what it's supposed to do.

By the way, I have not been able to find the matlab function tridieig. It doesn't seem to be a standard matlab function. Do you know what that's all about?

from nitime.

it seems that you found the good references ! so

option 1 : wrap the relevant function in scipy
option 2 : cythonize the gist

in recent matlab the function is tridisolve in the signal processing toolbox.

from nitime.

Hmm - those options sound hard. I think that the better option is probably #1.

As for understanding how Matlab does this, I haven't been able to find any functions in the signal processing toolbox starting with trid*

http://www.mathworks.com/help/toolbox/signal/f9-131178c.html

Do you think this is the one?

http://www.mathworks.com/matlabcentral/fileexchange/4822-using-numerical-computing-with-matlab-in-the-classroom/content/tridisolve.m

Doesn't seem like there's much there

from nitime.

the pb with #1 is that it requires me to work with the dev version of scipy which is a pain as I use EPD compiled with MKL for speed.

regarding the matlab function it is a private mex functions that's why you don't see it right away.

looks like we need a serious sprint on nitime :)

from nitime.

I don't have the matlab dpss code in front of me, but I think the answer may be interpolation. I think all the underlying wavefunctions are the same at any discretization, so you can compute a more course dpss and interpolate to the higher res one.

from nitime.

Mike - I am not sure that's the way matlab solves this problem, but it
seems like a very clever solution and might just work (and we could
test the solution against a matlab-generated solution). Alex - what do
you think? Do you see any problem with this solution? I am sorry for
not responding to this thread earlier. I have been swamped and haven't
really had time to work on nitime in the past few weeks. Does any one
of you have time to take a crack at it? Otherwise, I might be able to
make some time on Monday to try implementing this.

On Fri, May 20, 2011 at 12:31 PM, miketrumpis
[email protected]
wrote:

I don't have the matlab dpss code in front of me, but I think the answer may be interpolation. I think all the underlying wavefunctions are the same at any discretization, so you can compute a more course dpss and interpolate to the higher res one.

Reply to this email directly or view it on GitHub:
#61 (comment)

from nitime.

I really doubts Matlab uses interpolation. However this hack might work even if it's hack, hence it is not very satisfying. Can cython and C compiler be an acceptable dependency for nitime? If so I think we should consider cythonizing the code found by ariel.

from nitime.

OK - I have looked around a bit more. It turns out that matlab does use interpolation, at least in some cases. Their documentation also mentions that interpolation causes errors on the order of 10e-5 relative to direct calculation. I am going to start by implementing that and see how it works. I think that we can take on cython as a dependency, though I am a bit wary of having the kind of issues that I have seen dipy have, where things may or may not depending on the exact version of cython that is installed. It does open the possibility of using cython to optimize other parts of the code that we have already been considering cythonizing. So yes - if you want to work on a cython solution to this problem in parallel to my working on the interpolation, that would be great.

from nitime.

I don't want to push for cython (I try to avoid it in MNE) but most of the deployment pbs disappear if you accept to ship the generated C code and not just the .pyx file. I honestly do not have much time these days but I'll continue to review your patches and give you feedback.

from nitime.

OK - here's that solution, for now.

Alex - I am also a bit swamped these days, but I am trying to get to all your comments and suggestions. Your feedback is much appreciated.

from nitime.

I also addressed this issue in my last commit--sorry I didn't know how to link up the issue. See the dpss stuff in the commit

miketrumpis@6afdc1b

from nitime.

Following up on this, these commits integrate both solutions (interpolation and direct calculation) + some tests to make sure that both methods give the same result and that both are the same as a matlab-generated long window:

arokem@8819ea1

We'll merge it all into the trunk later today and I think we can close this one.

from nitime.