pranaysy / riksneurotools

Rik Henson's neuroimaging analysis tools
GNU General Public License v3.0
1 stars 0 forks source link

Norm used for scaling CL is neither Frobenius nor L2 #1

Open pranaysy opened 1 year ago

pranaysy commented 1 year ago

Problem

In fuse_lfp:L137 the $\text{norm}$ of the leadfield that has been reduced from $sensors\times vertices$ to $sensors\times ROIs$ does not correspond to any known $\text{norm}$.

More precisely, the $\text{norm}$ is estimated in code as: sqrt(trace(CL{m}*CL{m}')/Nchn(m)),
but, for matrices $\text{Frobenius norm}$ is defined as `sqrt(trace(CL{m}CL{m}'))(of many sources, [here's one at MATLAB's docs](https://uk.mathworks.com/help/matlab/ref/norm.html#bvhji30-5)) while $L_2 \text{ norm}$ for a matrix (also called the ***spectral norm***) is defined asmax(svd(CL{m}))` (MATLAB's docs)


For matrices: $\text{Frobenius norm}$ = norm(CL{m}, 'fro') = `sqrt(trace(CL{m}CL{m}'))=sqrt(sum(CL{m}(:).^2))`

Reproducible example A

x = ones(5);

% Frobenius norms
norm(x, 'fro')      % returns 5
sqrt(trace(x*x'))   % returns 5
sqrt(sum(x(:).^2))  % returns 5

% L2 norms
norm(x, 2)          % returns 5
max(svd(x))         % returns 5

% Implemented norm
sqrt(trace(x*x')/5) % returns 2.2361

Reproducible example B

x = magic(5);

% x =
%
%    17    24     1     8    15
%    23     5     7    14    16
%     4     6    13    20    22
%    10    12    19    21     3
%    11    18    25     2     9

% Frobenius norms
norm(x, 'fro')      % returns 74.330
sqrt(trace(x*x'))   % returns 74.330
sqrt(sum(x(:).^2))  % returns 74.330

% L2 norms
norm(x, 2)          % returns 65
max(svd(x))         % returns 65

% Implemented norm
sqrt(trace(x*x')/5) % returns 33.242

Solution

As the implemented $\text{norm}$ closely matches the $\text{Frobenius norm}$, switch to MATLAB's default implementation with norm(x, 'fro')

pranaysy commented 1 year ago

Lines 77-80 of spm_eeg_invert.m describe how scaling is performed for both the lead field and data.

% Potential scaling differences between the lead-fields are handled by
% scaling L{1} such that trace(L{1}*L{1}') = constant (number of spatial
% modes or channels), while scaling the data such that trace(AY{n}*AY{n}') =
% constant over subjects (and modalities; see below).

L316 of spm_eeg_invert.m performs this normalization, which is exactly the same as fuse_lfp.m#L137

% normalise lead-field
%------------------------------------------------------------------
Scale = sqrt(trace(UL{m}*UL{m}')/Nm(m));
UL{m} = UL{m}/Scale;
pranaysy commented 1 year ago

The normalization issue is addressed in this paper:

Henson, Richard N., Elias Mouchlianitis, and Karl J. Friston. “MEG and EEG Data Fusion: Simultaneous Localisation of Face-Evoked Responses.” NeuroImage 47, no. 2 (August 15, 2009): 581–89. https://doi.org/10.1016/j.neuroimage.2009.04.063.

In section An empirical Bayesian approach to MEG/EEG fusion:

image

This essentially argues that the relative sensitivity of the lead field may not be identical between different modalities. One way to normalize this is to make sure that the variance/scaling done by each sensor through the lead field is normalized to one. The trace of the covariance matrix of the lead field, $tr(LL^T)$, estimates the total variance/scaling done when $L$ is applied. The square root of of this trace yields the total standard deviation instead of total variance. So by dividing the lead field by $\sqrt{tr(LL^T)}$ or the $\text{Frobenius norm}$ ensures that the scaling achieved by $L$ is normalized to the $\text{Frobenius norm}$ of $L$ which is equivalent to the $\text{RMS gain}$ of $L$ and of course need not be one. However when the trace of the covariance matrix itself is normalized by the number of elements in the trace, $m$ (or on the diagonal), i.e. $\frac{1}{m} tr(LL^T)$, the meaning of this quantity changes from total scaling to the average scaling achieved when $L$ is applied. When $L$ is normalized by the square root of this quantity, it implies that the average scaling done by each row of $L$ is one. If there are $m$ elements in $L$ then the total scaling will therefore be $m$. This normalization procedure therefore acts in such a way that the average scaling obtained by each sensor regardless of modality will be one (exactly one if no sensor noise and iid sources)

Attached PDF: MEG and EEG data fusion - Simultaneous localisation of face-evoked responses - Henson et al [2009][NeuroImage].pdf


External links for further reading: