mne-python

Downloads: 1

[074d3d]: / doc / _includes / inverse.rst

History

Download this file

608 lines (446 with data), 26.2 kB

The minimum-norm current estimates

This section describes the mathematical details of the calculation of minimum-norm estimates. In Bayesian sense, the ensuing current distribution is the maximum a posteriori (MAP) estimate under the following assumptions:

The viable locations of the currents are constrained to the cortex. Optionally, the current orientations can be fixed to be normal to the cortical mantle.
The amplitudes of the currents have a Gaussian prior distribution with a known source covariance matrix.
The measured data contain additive noise with a Gaussian distribution with a known covariance matrix. The noise is not correlated over time.

Computing the inverse operator is accomplished using :func:`mne.minimum_norm.make_inverse_operator` and :func:`mne.minimum_norm.apply_inverse`. The use of these functions is presented in the tutorial :ref:`tut-inverse-methods`.

The linear inverse operator

The measured data in the source estimation procedure consists of MEG and EEG data, recorded on a total of N channels. The task is to estimate a total of Q strengths of sources located on the cortical mantle. If the number of source locations is P, Q = P for fixed-orientation sources and Q = 3P if the source orientations are unconstrained. The regularized linear inverse operator following from regularized maximal likelihood of the above probabilistic model is given by the Q×N matrix

M = R’G^⊤(GR’G^⊤ + C)^− 1 ,

where G is the gain matrix relating the source strengths to the measured MEG/EEG data, C is the data noise-covariance matrix and R’ is the source covariance matrix. The dimensions of these matrices are N×Q, N×N, and Q×Q, respectively. The Q×1 source-strength vector is obtained by multiplying the Q×1 data vector by Q.

The expected value of the current amplitudes at time t is then given by ĵ(t) = Mx(t), where x(t) is a vector containing the measured MEG and EEG data values at time t.

For computational convenience, the linear inverse operator is not computed explicitly. See :ref:`mne_solution` for mathematical details, and :ref:`CIHCFJEI` for a detailed example.

Regularization

The a priori variance of the currents is, in practice, unknown. We can express this by writing R’ = R ⁄ λ² = Rλ^− 2, which yields the inverse operator

M = R’G^⊤(GR’G^⊤ + C)^− 1 = Rλ^− 2G^⊤(GRλ^− 2G^⊤ + C)^− 1 = Rλ^− 2G^⊤λ²(GRG^⊤ + λ²C)^− 1 = RG^⊤(GRG^⊤ + λ²C)^− 1 ,

where the unknown current amplitude is now interpreted in terms of the regularization parameter λ². Larger λ² values correspond to spatially smoother and weaker current amplitudes, whereas smaller λ² values lead to the opposite.

We can arrive at the regularized linear inverse operator also by minimizing a cost function S with respect to the estimated current ĵ (given the measurement vector x at any given time t) as

min_ĵ{S} = min_ĵ{ẽ^⊤ẽ + λ²ĵ^⊤R^− 1ĵ} = min_ĵ{(x − Gĵ)^⊤C^− 1(x − Gĵ) + λ²ĵ^⊤R^− 1ĵ}

where the first term consists of the difference between the whitened measured data (see :ref:`whitening_and_scaling`) and those predicted by the model while the second term is a weighted-norm of the current estimate. It is seen that, with increasing λ², the source term receive more weight and larger discrepancy between the measured and predicted data is tolerable.

Whitening and scaling

The MNE software employs data whitening so that a 'whitened' inverse operator assumes the form

M̃ = MC^{¹ ⁄ ₂} = RG̃^⊤(G̃RG̃^⊤ + λ²I)^− 1 ,

where

G̃ = C^{− ¹ ⁄ ₂}G

is the spatially whitened gain matrix. We arrive at the whitened inverse operator equation :eq:`inv_m_tilde` by making the substitution for G from :eq:`inv_g_tilde` in :eq:`inv_m` as

M̃ = MC^{¹ ⁄ ₂} = RG^⊤(GRG^⊤ + λ²C)^− 1C^{¹ ⁄ ₂} = RG̃^⊤C^{¹ ⁄ ₂}(C^{¹ ⁄ ₂}G̃RG̃^⊤C^{¹ ⁄ ₂} + λ²C)^− 1C^{¹ ⁄ ₂} = RG̃^⊤C^{¹ ⁄ ₂}(C^{¹ ⁄ ₂}(G̃RG̃^⊤ + λ²I)C^{¹ ⁄ ₂})^− 1C^{¹ ⁄ ₂} = RG̃^⊤C^{¹ ⁄ ₂}C^{− ¹ ⁄ ₂}(G̃RG̃^⊤ + λ²I)^− 1C^{− ¹ ⁄ ₂}C^{¹ ⁄ ₂} = RG̃^⊤(G̃RG̃^⊤ + λ²I)^− 1 .

The expected current values are

ĵ(t) = Mx(t) = MC^{¹ ⁄ ₂}C^{− ¹ ⁄ ₂}x(t) = M̃x̃(t)

knowing :eq:`inv_m_tilde` and taking

x̃(t) = C^{− ¹ ⁄ ₂}x(t)

as the whitened measurement vector at time t. The spatial whitening operator C^{− ¹ ⁄ ₂} is obtained with the help of the eigenvalue decomposition C = U_CΛ²_CU^⊤_C as C^{− ¹ ⁄ ₂} = Λ^− 1_CU^⊤_C. In the MNE software the noise-covariance matrix is stored as the one applying to raw data. To reflect the decrease of noise due to averaging, this matrix, C₀, is scaled by the number of averages, L, i.e., C = C₀ ⁄ L.

Note

When EEG data are included, the gain matrix G needs to be average referenced when computing the linear inverse operator M. This is incorporated during creating the spatial whitening operator C^{− ¹ ⁄ ₂}, which includes any projectors on the data. EEG data average reference (using a projector) is mandatory for source modeling and is checked when calculating the inverse operator.

As shown above, regularization of the inverse solution is equivalent to a change in the variance of the current amplitudes in the Bayesian a priori distribution.

A convenient choice for the source-covariance matrix R is such that trace(G̃RG̃^⊤) ⁄ trace(I) = 1. With this choice we can approximate λ² ∼ 1 ⁄ \rmSNR², where SNR is the (amplitude) signal-to-noise ratio of the whitened data.

Note

The definition of the signal to noise-ratio/ λ² relationship given above works nicely for the whitened forward solution. In the un-whitened case scaling with the trace ratio trace(GRG^⊤) ⁄ trace(C) does not make sense, since the diagonal elements summed have, in general, different units of measure. For example, the MEG data are expressed in T or T/m whereas the unit of EEG is Volts.

See :ref:`tut-compute-covariance` for example of noise covariance computation and whitening.

Regularization of the noise-covariance matrix

Since finite amount of data is usually available to compute an estimate of the noise-covariance matrix C, the smallest eigenvalues of its estimate are usually inaccurate and smaller than the true eigenvalues. Depending on the seriousness of this problem, the following quantities can be affected:

The model data predicted by the current estimate,
Estimates of signal-to-noise ratios, which lead to estimates of the required regularization, see :ref:`mne_regularization`,
The estimated current values, and
The noise-normalized estimates, see :ref:`noise_normalization`.

Fortunately, the latter two are least likely to be affected due to regularization of the estimates. However, in some cases especially the EEG part of the noise-covariance matrix estimate can be deficient, i.e., it may possess very small eigenvalues and thus regularization of the noise-covariance matrix is advisable.

Historically, the MNE software accomplishes the regularization by replacing a noise-covariance matrix estimate C with

C’ = C + ∑_kε_k^‒σ_k²I^(k) ,

where the index k goes across the different channel groups (MEG planar gradiometers, MEG axial gradiometers and magnetometers, and EEG), ε_k are the corresponding regularization factors, ^‒σ_k are the average variances across the channel groups, and I^(k) are diagonal matrices containing ones at the positions corresponding to the channels contained in each channel group.

See :ref:`plot_compute_covariance_howto` for details on computing and regularizing the channel covariance matrix.

Computation of the solution

The most straightforward approach to calculate the MNE is to employ the expression of the original or whitened inverse operator directly. However, for computational convenience we prefer to take another route, which employs the singular-value decomposition (SVD) of the matrix

A = G̃R^{¹ ⁄ ₂} = UΛV^⊤

where the superscript ¹ ⁄ ₂ indicates a square root of R. For a diagonal matrix, one simply takes the square root of R while in the more general case one can use the Cholesky factorization R = R_CR^⊤_C and thus R^{¹ ⁄ ₂} = R_C.

Combining the SVD from :eq:`inv_a` with the inverse equation :eq:`inv_m` it is easy to show that

M̃ = RG̃^⊤(G̃RG̃^⊤ + λ²I)^− 1 = R^{¹ ⁄ ₂}A^⊤(AA^⊤ + λ²I)^− 1 = R^{¹ ⁄ ₂}VΛU^⊤(UΛV^⊤VΛU^⊤ + λ²I)^− 1 = R^{¹ ⁄ ₂}VΛU^⊤(U(Λ² + λ²I)U^⊤)^− 1 = R^{¹ ⁄ ₂}VΛU^⊤U(Λ² + λ²I)^− 1U^⊤ = R^{¹ ⁄ ₂}VΛ(Λ² + λ²I)^− 1U^⊤ = R^{¹ ⁄ ₂}VΓU^⊤

where the elements of the diagonal matrix Γ are simply

γ_k = (λ_k)/(λ²_k + λ²) .

From our expected current equation :eq:`inv_j_hat_t` and our whitened measurement equation :eq:`inv_tilde_x_t`, if we take

w(t) = U^⊤x̃(t) = U^⊤C^{− ¹ ⁄ ₂}x(t) ,

we can see that the expression for the expected current is just

ĵ(t) = R^{¹ ⁄ ₂}VΓw(t) = ∑_k^‒v_kγ_kw_k(t) ,

where ^‒v_k = R^{¹ ⁄ ₂}v_k, with v_k being the k th column of V. It is thus seen that the current estimate is a weighted sum of the "weighted" eigenleads v_k.

It is easy to see that w(t) ∝ √(L). To maintain the relation (G̃RG̃^⊤) ⁄ trace(I) = 1 when L changes we must have R ∝ 1 ⁄ L. With this approach, λ_k is independent of L and, for fixed λ, we see directly that j(t) is independent of L.

The minimum-norm estimate is computed using this procedure in :func:`mne.minimum_norm.make_inverse_operator`, and its usage is illustrated in :ref:`CIHCFJEI`.

Noise normalization

Noise normalization serves three purposes:

It converts the expected current value into a dimensionless statistical test variable. Thus the resulting time and location dependent values are often referred to as dynamic statistical parameter maps (dSPM).
It reduces the location bias of the estimates. In particular, the tendency of the MNE to prefer superficial currents is eliminated.
The width of the point-spread function becomes less dependent on the source location on the cortical mantle. The point-spread is defined as the MNE resulting from the signals coming from a point current source (a current dipole) located at a certain point on the cortex.

In practice, noise normalization is implemented as a division by the square root of the estimated variance of each voxel. In computing these noise normalization factors, it's convenient to reuse our "weighted eigenleads" definition from equation :eq:`inv_j_hat_t` in matrix form as

^‒V = R^{¹ ⁄ ₂}V .

dSPM

Noise-normalized linear estimates introduced by Dale et al. :footcite:`DaleEtAl1999` require division of the expected current amplitude by its variance. In practice, this requires the computation of the diagonal elements of the following matrix, using SVD equation :eq:`inv_m_tilde` and :eq:`inv_eigenleads_weighted`:

MCM^⊤ = MC^{¹ ⁄ ₂}C^{¹ ⁄ ₂}M^⊤ = M̃M̃^⊤ = R^{¹ ⁄ ₂}VΓU^⊤UΓV^⊤R^{¹ ⁄ ₂} = ^‒VΓ²^‒V^⊤ .

Because we only care about the diagonal entries here, we can find the variances for each source as

σ²_k = γ²_k

Under the conditions expressed at the end of :ref:`mne_solution`, it follows that the t-statistic values associated with fixed-orientation sources) are thus proportional to √(L) while the F-statistic employed with free-orientation sources is proportional to L, correspondingly.

Note

The MNE software usually computes the square roots of the F-statistic to be displayed on the inflated cortical surfaces. These are also proportional to √(L).

sLORETA

sLORETA :footcite:`Pascual-Marqui2002` estimates the current variances as the diagonal entries of the resolution matrix, which is the product of the inverse and forward operators. In other words, the diagonal entries of (using :eq:`inv_m_tilde_svd`, :eq:`inv_g_tilde`, and :eq:`inv_a`)

MG = MC^{¹ ⁄ ₂}C^{− ¹ ⁄ ₂}G = M̃G̃ = R^{¹ ⁄ ₂}VΓU^⊤G̃R^{¹ ⁄ ₂}R^{− ¹ ⁄ ₂} = R^{¹ ⁄ ₂}VΓU^⊤UΛV^⊤R^{− ¹ ⁄ ₂} = R^{¹ ⁄ ₂}VΓU^⊤UΛV^⊤R^{¹ ⁄ ₂}R^− 1 = ^‒VΓU^⊤UΛ^‒V^⊤R^− 1 = ^‒VΓΛ^‒V^⊤R^− 1 .

Because R is diagonal and we only care about the diagonal entries, we can find our variance estimates as

σ²_k = γ_kλ_kR^− 1_k, k = ⎛⎝(λ_k)/((λ²_k + λ²))⎞⎠⎛⎝(λ_k)/(1)⎞⎠⎛⎝(1)/(λ²)⎞⎠ = (λ²_k)/((λ²_k + λ²)λ²) = ⎛⎝(λ²_k)/((λ²_k + λ²)²)⎞⎠⎛⎝(λ² + λ²_k)/(λ²)⎞⎠ = ⎛⎝(λ_k)/(λ²_k + λ²)⎞⎠²⎛⎝1 + (λ²_k)/(λ²)⎞⎠ = γ²_k⎛⎝1 + (λ²_k)/(λ²)⎞⎠ .

eLORETA

While dSPM and sLORETA solve for noise normalization weights σ²_k that are applied to standard minimum-norm estimates ĵ(t), eLORETA :footcite:`Pascual-Marqui2011` instead solves for a source covariance matrix R that achieves zero localization bias. For fixed-orientation solutions the resulting matrix R will be a diagonal matrix, and for free-orientation solutions it will be a block-diagonal matrix with 3×3 blocks.

In :footcite:`Pascual-Marqui2011` eq. 2.13 states that the following system of equations can be used to find the weights, ∀i ∈ 1, ..., P (note that here we represent the equations from that paper using our notation):

r_i = [G^⊤_i(GRG^⊤ + λ²C)^− 1G_i]^{− ¹ ⁄ ₂}

And an iterative algorithm can be used to find the values for the weights r_i that satisfy these equations as:

Initialize identity weights.
Compute N = (GRG^⊤ + λ²C)^− 1.
Holding N fixed, compute new weights r_i = [G^⊤_iNG_i]^{− ¹ ⁄ ₂}.
Using new weights, go to step (2) until convergence.

In particular, for step (2) we can use our substitution from :eq:`inv_g_tilde` as:

N = (GRG^⊤ + λ²C)^− 1 = (C^{¹ ⁄ ₂}G̃RG̃^⊤C^{¹ ⁄ ₂} + λ²C)^− 1 = (C^{¹ ⁄ ₂}(G̃RG̃^⊤ + λ²I)C^{¹ ⁄ ₂})^− 1 = C^{− ¹ ⁄ ₂}(G̃RG̃^⊤ + λ²I)^− 1C^{− ¹ ⁄ ₂} = C^{− ¹ ⁄ ₂}(G̃RG̃^⊤ + λ²I)^− 1C^{− ¹ ⁄ ₂} .

Then defining Ñ as the whitened version of N, i.e., the regularized pseudoinverse of G̃RG̃^⊤, we can compute N as:

N = C^{− ¹ ⁄ ₂}(U_G̃RG̃^⊤Λ_G̃RG̃^⊤V^⊤_G̃RG̃^⊤ + λ²I)^− 1C^{− ¹ ⁄ ₂} = C^{− ¹ ⁄ ₂}(U_G̃RG̃^⊤(Λ_G̃RG̃^⊤ + λ²I)V^⊤_G̃RG̃^⊤)^− 1C^{− ¹ ⁄ ₂} = C^{− ¹ ⁄ ₂}V_G̃RG̃^⊤(Λ_G̃RG̃^⊤ + λ²I)^− 1U^⊤_G̃RG̃^⊤C^{− ¹ ⁄ ₂} = C^{− ¹ ⁄ ₂}ÑC^{− ¹ ⁄ ₂} .

In step (3) we left and right multiply with subsets of G, but making the substitution :eq:`inv_g_tilde` we see that we equivalently compute:

r_i = [G^⊤_iNG_i]^{− ¹ ⁄ ₂} = [(C^{¹ ⁄ ₂}G̃_i)^⊤NC^{¹ ⁄ ₂}G̃_i]^{− ¹ ⁄ ₂} = [G̃^⊤_iC^{¹ ⁄ ₂}NC^{¹ ⁄ ₂}G̃_i]^{− ¹ ⁄ ₂} = [G̃^⊤_iC^{¹ ⁄ ₂}C^{− ¹ ⁄ ₂}ÑC^{− ¹ ⁄ ₂}C^{¹ ⁄ ₂}G̃_i]^{− ¹ ⁄ ₂} = [G̃^⊤_iÑG̃_i]^{− ¹ ⁄ ₂} .

For convenience, we thus never need to compute N itself but can instead compute the whitened version Ñ.

Predicted data

Under noiseless conditions the SNR is infinite and thus leads to λ² = 0 and the minimum-norm estimate explains the measured data perfectly. Under realistic conditions, however, λ² > 0 and there is a misfit between measured data and those predicted by the MNE. Comparison of the predicted data, here denoted by x(t), and measured one can give valuable insight on the correctness of the regularization applied.

In the SVD approach we easily find

x̂(t) = Gĵ(t) = C^{¹ ⁄ ₂}UΠw(t) ,

where the diagonal matrix Π has elements π_k = λ_kγ_k The predicted data is thus expressed as the weighted sum of the 'recolored eigenfields' in C^{¹ ⁄ ₂}U.

Cortical patch statistics

If the add_dists=True option was used in source space creation, the source space file will contain Cortical Patch Statistics (CPS) for each vertex of the cortical surface. The CPS provide information about the source space point closest to it as well as the distance from the vertex to this source space point. The vertices for which a given source space point is the nearest one define the cortical patch associated with with the source space point. Once these data are available, it is straightforward to compute the following cortical patch statistics for each source location d:

The average over the normals of at the vertices in a patch, ^‒n_d,
The areas of the patches, A_d, and
The average deviation of the vertex normals in a patch from their average, σ_d, given in degrees.

use_cps parameter in :func:`mne.convert_forward_solution`, and :func:`mne.minimum_norm.make_inverse_operator` controls whether to use cortical patch statistics (CPS) to define normal orientations or not (see :ref:`CHDBBCEJ`).

Orientation constraints

The principal sources of MEG and EEG signals are generally believed to be postsynaptic currents in the cortical pyramidal neurons. Since the net primary current associated with these microscopic events is oriented normal to the cortical mantle, it is reasonable to use the cortical normal orientation as a constraint in source estimation. In addition to allowing completely free source orientations, the MNE software implements three orientation constraints based of the surface normal data:

Source orientation can be rigidly fixed to the surface normal direction by specifying fixed=True in :func:`mne.minimum_norm.make_inverse_operator`. If cortical patch statistics are available the average normal over each patch, ^‒n_d, are used to define the source orientation. Otherwise, the vertex normal at the source space location is employed.
A location independent or fixed loose orientation constraint (fLOC) can be employed by specifying fixed=False and loose=1.0 when calling :func:`mne.minimum_norm.make_inverse_operator` (see :ref:`plot_dipole_orientations_fLOC_orientations`). In this approach, a source coordinate system based on the local surface orientation at the source location is employed. By default, the three columns of the gain matrix G, associated with a given source location, are the fields of unit dipoles pointing to the directions of the x, y, and z axis of the coordinate system employed in the forward calculation (usually the :ref:`MEG head coordinate frame <head_device_coords>`). For LOC the orientation is changed so that the first two source components lie in the plane normal to the surface normal at the source location and the third component is aligned with it. Thereafter, the variance of the source components tangential to the cortical surface are reduced by a factor defined by the --loose option.
A variable loose orientation constraint (vLOC) can be employed by specifying fixed=False and loose parameters when calling :func:`mne.minimum_norm.make_inverse_operator` (see :ref:`plot_dipole_orientations_vLOC_orientations`). This is similar to fLOC except that the value given with the loose parameter will be multiplied by σ_d, defined above.

Depth weighting

The minimum-norm estimates have a bias towards superficial currents. This tendency can be alleviated by adjusting the source covariance matrix R to favor deeper source locations. In the depth weighting scheme employed in MNE analyze, the elements of R corresponding to the p th source location are be scaled by a factor

f_p = (g^⊤_1pg_1p + g^⊤_2pg_2p + g^⊤_3pg_3p)^− γ ,

where g_1p, g_2p, and g_3p are the three columns of G corresponding to source location p and γ is the order of the depth weighting, which is specified via the depth option in :func:`mne.minimum_norm.make_inverse_operator`.

Effective number of averages

It is often the case that the epoch to be analyzed is a linear combination over conditions rather than one of the original averages computed. As stated above, the noise-covariance matrix computed is originally one corresponding to raw data. Therefore, it has to be scaled correctly to correspond to the actual or effective number of epochs in the condition to be analyzed. In general, we have

C = C₀ ⁄ L_eff

where L_eff is the effective number of averages. To calculate L_eff for an arbitrary linear combination of conditions

y(t) = ⁿ∑_i = 1w_ix_i(t)

we make use of the the fact that the noise-covariance matrix

C_y = ⁿ∑_i = 1w²_iC_{x_i} = C₀ⁿ∑_i = 1w²_i ⁄ L_i

which leads to

1 ⁄ L_eff = ⁿ∑_i = 1w²_i ⁄ L_i

An important special case of the above is a weighted average, where

w_i = L_i ⁄ ⁿ∑_i = 1L_i

and, therefore

L_eff = ⁿ∑_i = 1L_i

Instead of a weighted average, one often computes a weighted sum, a simplest case being a difference or sum of two categories. For a difference w₁ = 1 and w₂ = − 1 and thus

1 ⁄ L_eff = 1 ⁄ L₁ + 1 ⁄ L₂

L_eff = (L₁L₂)/(L₁ + L₂)

Interestingly, the same holds for a sum, where w₁ = w₂ = 1. Generalizing, for any combination of sums and differences, where w_i = 1 or w_i = − 1, i = 1…n, we have

1 ⁄ L_eff = ⁿ∑_i = 11 ⁄ L_i