Hi Jeff,

You are trying to use a classic “imaging resolution” approach to this kind of data, which in general does not work. Pointing at one of our papers,

S. Baillet, JC Mosher, RM Leahy, Electromagnetic Brain Mapping, IEEE Signal Processing Magazine 18 (6): 14-30, pp 18 Nov 2001. http://dx.doi.org/10.1109/79.962275 (see also http://neuroimage.usc.edu/brainstorm/Pub)

Equation (23) as an example. Let’s eliminate considerations of “W” by setting it to the identity matrix, which as you point out is just one of the many flavors of minimum norms (depth weighting, LORETA, LAURA, etc), thus we are estimating [B]S[/B], the time series, to which you want to reapply [B]A[/B], the lead field matrix, and reconstruct the original data.

In your example, [B]A[/B] is the 256 x 2403 lead field matrix, for 256 channels of EEG data (256 lead fields) evaluated at 2,403 fixed orientation dipolar points in the brain volume. In Equation (23), we have set the noise covariance matrix to identity , with the assumption that the user has already pre-whitened the modeling equations, thus the remaining “noise” parameter is “lambda”, also known as the Tikhonov regularization parameter.

In noiseless data (which does not exist, but let’s pretend), then we set lambda to zero, and our estimate of the 2,403 time series is then

[B]S[/B][SUB]e[/SUB][SUP]T[/SUP] = [B]A[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP])[SUP]-1[/SUP][B]M[/B]

where [B]M[/B] are your noiseless measurements. Thus, in a perfect world,

[B]D[/B][SUB]e [/SUB] = [B]A[/B] [B]S[/B][SUB]e[/SUB][SUP]T[/SUP] = [B]AA[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP])[SUP]-1[/SUP][B]M[/B]

and thus in a perfect world, then [B]D[/B][SUB]e [/SUB] = [B]M[/B], which is what you are expecting.

Note, however, that (1) it would not matter at all what [B]A[/B] you used, you would always get a perfect estimate, and therefore we have not truly tested the correctness of the calculations in [B]A[/B], which was your original goal. This fault in testing is known as the “inverse crime” where you use the same model matrix in both the forward and inverse procedures. You could just as easily create random matrices of [B]A[/B]. To avoid the inverse crime, you need an independently derived forward model [B]A[/B][SUB]2[/SUB] for the reconstruction,

[B]D[/B][SUB]e [/SUB] = [B]A[/B][SUB]2 [/SUB] [B]S[/B][SUB]e[/SUB][SUP]T[/SUP] = [B]A[SUB]2[/SUB]A[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP])[SUP]-1[/SUP][B]M[/B].

(2) This reconstruction requires the inversion ([B]AA[/B][SUP]T[/SUP])[SUP]-1[/SUP] in order to form the initial time series estimate. Look at the eigenvalues of this “Gramian” [B]AA[/B][SUP]T[/SUP] in your lead field model, of the 256 eigenvalues, the last several will be extremely small, many orders of magnitude smaller than is physically plausible, and thus generally rejected in the estimation process. We control this range of rejections through the parameters of noise covariance, whitening, Signal to Noise Ratio. If you set our SNR parameter to 1e15, for example, turn off whitening, turn off regularization, you should find most of your sources passing undisturbed through the estimation, but we may still eliminate some leadfields simply because we still detect them as unstable. But note (1) above, this is overall not a good test anyway

In our equation (23), we control the instability of the Gramian through the parameter “lambda”, which in a Bayesian framework is tied to the SNR, or may simply be considered a Tikhonov regularizer in a more general framework. It automatically suppresses the tiny eigenvalues, so that the estimate is

[B]S[/B][SUB]e[/SUB][SUP]T[/SUP] = [B]A[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP] + “lam” [B]I[/B])[SUP]-1[/SUP][B]M[/B]

and thus

[B]D[/B][SUB]e [/SUB] = [B]A[/B] [B]S[/B][SUB]e[/SUB][SUP]T[/SUP] = [B]AA[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP]+ “lam” [B]I[/B])[SUP]-1[/SUP][B]M[/B]

The “imaging resolution kernel” is now [B]AA[/B][SUP]T[/SUP] ([B]AA[/B][SUP]T[/SUP]+ “lam” [B]I[/B])[SUP]-1[/SUP], which in general cannot reconstruct all EEG channels.

The theoretical issue is that the Gramian in quasistatic electromagnetics is fundamentally limited in the number of significant eigenvalues, i.e you will never be able to reconstruct all 256 lead fields in your case in any practical sense.

– John