Back to the main page.
Bug 2838 - Sourceanalysis using MNE doesn't use noise covariance
| Status | CLOSED FIXED |
| Reported | 2015-02-11 14:53:00 +0100 |
| Modified | 2015-07-15 13:31:25 +0200 |
| Product: | FieldTrip |
| Component: | inverse |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Linux |
| Importance: | P5 normal |
| Assigned to: | Jan-Mathijs Schoffelen |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Pim Mostert - 2015-02-11 14:53:36 +0100
cfg = []; cfg.method = 'mne'; cfg.grid = lf; cfg.grid.pos = sourcespace.pnt; cfg.vol = vol; cfg.mne.lambda = 5; cfg.mne.keepfilter = 'yes'; sourceDiff = ft_sourceanalysis(cfg, diff); Where diff was obtained using ft_timelockanalysis and cfg.covariance = 'yes'. The verbal output reports that the noise covariance is used. However, I tried changing the covariance matrix of diff in many different ways [even diff.cov = rand(size(diff.cov))], but the results are always the same. Also, changing cfg.mne.lambda doesn't make a change. I took the liberty to look into the code, and it appears that in minimumnormestimate.m, lines 242-246, the noise covariance matrix indeed isn't used at all. Thanks, Pim
Jan-Mathijs Schoffelen - 2015-02-11 16:00:54 +0100
You're absolutely right, I think! It looks as if the code is only correct when you specify cfg.mne.prewhiten = 'yes', in which case the leadfields are prewhitened with C. Note that the prewhitening operation essentially reduces the prewhitened noise covariance to an identity matrix, in which case the lambda parameter can assume a value irrespective of the units that are in the leadfields and the data. (In other words, when you don't prewhiten (and provided the commented out code (line 232-238) is used) the relative scaling of C versus A*R*A' determines the magnitude of meaningful lambda values. I will adjust the code to run correctly when prewhiten is 'no'. Thanks for noticing!
Jan-Mathijs Schoffelen - 2015-02-11 16:13:48 +0100
[jansch@mentat003 inverse]$ svn diff minimumnormestimate.m Index: minimumnormestimate.m =================================================================== --- minimumnormestimate.m (revision 10204) +++ minimumnormestimate.m (working copy) @@ -228,26 +228,27 @@ lambda = trace(A * R * A')/(trace(C)*snr^2); end - %% equation 5 from Lin et al 2004 (this implements Dale et al 2000, and Liu et al. 2002) - %denom = (A*R*A'+(lambda^2)*C); - %if cond(denom)<1e12 - % w = R * A' / denom; - %else - % fprintf('taking pseudo-inverse due to large condition number\n'); - % w = R * A' * pinv(denom); - %end - - % as documented on MNE website, this is replacing the part of the code above, it gives - % more stable results numerically. - Rc = chol(R, 'lower'); - [U,S,V] = svd(A * Rc, 'econ'); - s = diag(S); - ss = s ./ (s.^2 + lambda); - w = Rc * V * diag(ss) * U'; - - % unwhiten the filters to bring them back into signal subspace - if dowhiten + if dowhiten, + % as documented on MNE website, this is replacing the part of the code above, it gives + % more stable results numerically. + Rc = chol(R, 'lower'); + [U,S,V] = svd(A * Rc, 'econ'); + s = diag(S); + ss = s ./ (s.^2 + lambda); + w = Rc * V * diag(ss) * U'; + + % unwhiten the filters to bring them back into signal subspace w = w*P; + + else + %% equation 5 from Lin et al 2004 (this implements Dale et al 2000, and Liu et al. 2002) + denom = (A*R*A'+(lambda^2)*C); + if cond(denom)<1e12 + w = R * A' / denom; + else + fprintf('taking pseudo-inverse due to large condition number\n'); + w = R * A' * pinv(denom); + end end end % if empty noisecov
Jan-Mathijs Schoffelen - 2015-02-11 16:14:34 +0100
[jansch@mentat003 inverse]$ svn commit -m "bugfix - properly invert the leadfield when prewhiten is set to 'no'" minimumnormestimate.m Sending minimumnormestimate.m Transmitting file data . Committed revision 10207.