% Here's my two cents ElisaM.
% Sorry this took so long. I'm just very busy, and I haven't been working
% on sLORETA for a while, even though I did write the sLORETA code for BST.
% I now use eLORETA (not sLORETA) for simple distributed linear imaging instead,
% and I'm working on COMPLETELY different methods for inverse imaging.
% This is an actual program written in Matlab tongue, so you can copy paste
% and reproduce everything as many times as you want.
L=randn(10,100); % Toy problem, "lead field" matrix is sampled from a Gaussian distribution, only 10 channels, 100 dipoles.
Ln=randn(10,13); % Add a noise dictionary to get covariance matrix with non-diagonal structure (i.e., correlatted noise).
Jn=randn(13,1000); % The noise current density time series associated with the noise dictionary.
Bn=LnJn; % Simulated noise data.
Cn=cov(Bn')/100; % Noise covariance matrix.
Cn=(Cn+Cn')/2; % Making sure it's perfectly symmetric.
W=inv(sqrtm(Cn)); % Whitening matrix.
LW=WL; % Whitened lead field.
Sigma_J=eye(100,100)(10/trace(LWLW')); % A priori source covariance.
% Numerically more stable SVD approach. I think this is how I wrote for
% Brainstorm many years ago.
Rc=chol(Sigma_J,'lower'); % Cholesky used only in case we had nonzero off-diagonal elements, which we usually never do, at least in BST.
LWRc=LWRc; % This is the matrix we need to SVD.
[U,S,V]=svd(LWRc,'econ'); % SVDing.
s=diag(S); % vector of singular values.
SNR=3; % Just a common safe default, could be time-dependent based on whitened data.
lambda2=SNR.^(-2); % Would be time dependent if SNR was time-dependent.
ss=s./(s.^2 + lambda2); % Same as above.
Kernel=RcV*diag(ss)U'; % This is the MNE "imaging kernel".
R=KernelLW; % Resolution matrix for MNE "imaging kernel"
sloretadiag=sqrt(sum(Kernel.*LW',2)); % computing the square root of diagonal of R above. Equivalently, sloretadiag=sqrt(diag(R));
KernelSLO=Kernel./repmat(sloretadiag(:),[1 10]); % sLORETA operator for the simpler case of fixed dipole orientations.
% Gram numerically noisy approach (i.e., inner product might accumulate error when there's a lot of dipoles, but it doesn't really matter in noisy reality)
T=Sigma_JLW'inv(LWSigma_JLW' + lamda2*eye(10,10)); % Exactly equivalent to MNE Kernel via SVD,
% which is the same as:
T_alt=(1/lambda2)Sigma_JLW'inv(LW(1/lambda2)Sigma_JLW' + eye(10,10)); % Exactly equivalent to MNE Kernel via SVD
corrcoef(Kernel,T)
norm(Kernel)
norm(T)
corrcoef(T,T_alt)
norm(T)
norm(T_alt)
T_oldschool=Sigma_JL'inv(LSigma_JL' + lamda2Cn); % Matrix is not identical to Kernel (as it should be) but it will produce the same
% exact current density vector when multiplied by data vector that you get if you multiply Kernel
% by whitened data vector.
% Note that L is used not LW.
jK=KernelLW(:,1); % essentialy the point-spread function for first dipole
j_oldschool=T_oldschool*L(:,1); % essentialy the point-spread function for first dipole
corrcoef(j,jK) % perfect correlation
norm(jK) % same norms
norm(j) % same norms
RT=TLW; % Matches R as it should. Resolution matrix obviously has to since T and Kernel are identical.
R_oldschool=T_oldschoolL; % Also matches R as it should. Note that you have to multiply by L not LW.
corrcoef(R,R_oldschool)
norm(R)
norm(R_oldschool)
% So all these different ways of computing the MNE resolution matrix are going
% to give you the same exact thing. Therefore, they will give you the same
% exact sLORETA 'inverse' operator. The important matrix here again, I
% want to emphasize is the resolution matrxix (as far as I know).
% Now let's look at equation (15) of Pascual-Marqui 2002
S_j_PM=T_oldschool*(LSigma_JL' + lambda2Cn)T_oldschool';
% Is it exactly equal to the resolution matrix? Let's check.
corrcoef(R,S_j_PM) % well you should get a perfect correlation of 1 (just a fast check as we hack along)
% But is there a scale difference?
norm(R)
norm(S_j_PM) % Yes it seems there's a scale factor that is off.
% To get the exact same scale for the resolution matrix.
% You have to do this:
S_j_alt=(1/Sigma_J(1,1))T_oldschool(LSigma_JL' + lambda2*Cn)*T_oldschool'; % Note that the factor is 1/Sigma_J(1,1) not lambda2 as in your email, these are different.
% Let's check correlation
corrcoef(R,S_j_alt)
norm(R)
norm(S_j_alt) % Bingo they match now, which is what I think you had noticed. You are very very veru observant.
% But what's important for sLORETA is the resolution matrix.
% Let's tale a look at the old school resolution matrix.
R_oldschool=T_oldschool*L; % Note that you multiply by L not LW in this case.
corrcoef(R,R_oldschool) % correlation of 1.
norm(R) % checking scales
norm(R_oldschool) % yes they match, and they are identical.
% So my two cents is that as far as I know what really matters for computing sLORETA
% is the MNE resolution matrix, which was done here in several different
% ways and gets the same exact result. You were right that there's a
% scaling issue with equation (15) but that doesn't matter as far as I
% know. You have to use the resolution matrix (see Backus-Gilbert work for
% more info on resoution matrices and their optimization).
% AS A FINAL NOTE, I REALLY WANT TO EMPHASIZE THAT I PERSONALLY STRONGLY
% RECOMMEND ELORETA OVER ANY OTHER LINEAR DISTRIBURED METHOD INCLUDING
% SLORETA. WITH ELORETA YOU GET ZERO LOCALIZATION BIAS AND THEY ARE
% CURRENT DENISTY SOLUTIONS, MEANING REAL FEASIBLE SOLUTIONS, NOT ACTIVITY
% MAPS WHICH IS WHAT SLORETA OR BEAMFORMERS GIVE YOU, YUCK.
% About question #2 I think you have to ask Francois but I think basically
% he probably multiplies the imaging Kernel by the whitener so that the
% data is used instead of the whitened data. If this is the case, it
% actually uses more memory and more floating point operations, but there
% might be a practical reason. And as for the EEG case, well I have a
% feeling that the data is already rereferenced to the average 'electrode' so no need to
% keep doing that. The data is already probably stored properly referenced to the mean.
% Hope it makes sense. Best regards. And keep looking for
% mistakes as there could be some waiting to be corrected.
% In any case, if there had been an error, it would have been a simple
% scaling and that wouldn't matter much for sLORETA as those activity maps
% are in arbitrary units, not units of current densisty.
% These are only my two cents. Anyone should feel more than free to look
% into this some more and possibly correct me if I'm wrong this time.
% It's very tricky business, so the more brains the merrier. Brainstorm
% away. –Rey Rene Ramirez