Skip to content

Linearly Constrained Minimum Variance

Solver ID: LCMV

Usage

from invert import Solver

# fwd = ...    (mne.Forward object)
# evoked = ... (mne.Evoked object)

solver = Solver("LCMV")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview

Classic time-domain linearly constrained minimum-variance (LCMV) beamformer / spatial filter.

References

  1. Van Veen, B. D., van Drongelen, W., Yuchtman, M., & Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Transactions on Biomedical Engineering, 44(9), 867-880.
  2. Van Veen, B. D., & Buckley, K. M. (1988). Beamforming: A versatile approach to spatial filtering. IEEE ASSP Magazine, 5(2), 4-24.

API Reference

Bases: BaseSolver

Class for the Linearly Constrained Minimum Variance Beamformer (LCMV) inverse solution.

References

[1] Van Veen, B. D., & Buckley, K. M. (1988). Beamforming: A versatile approach to spatial filtering. IEEE ASSP Magazine, 5(2), 4-24.

Source code in invert/solvers/beamformers/lcmv.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
class SolverLCMV(BaseSolver):
    """Class for the Linearly Constrained Minimum Variance Beamformer (LCMV) inverse solution.

    References
    ----------
    [1] Van Veen, B. D., & Buckley, K. M. (1988). Beamforming: A versatile
        approach to spatial filtering. IEEE ASSP Magazine, 5(2), 4-24.
    """

    meta = SolverMeta(
        slug="lcmv",
        full_name="Linearly Constrained Minimum Variance",
        category="Beamformers",
        description=(
            "Classic time-domain linearly constrained minimum-variance (LCMV) "
            "beamformer / spatial filter."
        ),
        references=[
            "Van Veen, B. D., van Drongelen, W., Yuchtman, M., & Suzuki, A. (1997). "
            "Localization of brain electrical activity via linearly constrained minimum "
            "variance spatial filtering. IEEE Transactions on Biomedical Engineering, "
            "44(9), 867-880.",
            "Van Veen, B. D., & Buckley, K. M. (1988). Beamforming: A versatile "
            "approach to spatial filtering. IEEE ASSP Magazine, 5(2), 4-24.",
        ],
    )

    def __init__(self, name="LCMV Beamformer", reduce_rank=True, rank="auto", **kwargs):
        self.name = name
        return super().__init__(reduce_rank=reduce_rank, rank=rank, **kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj,
        *args,
        alpha="auto",
        weight_norm=True,
        verbose=0,
        **kwargs,
    ):
        """Calculate inverse operator.

        Parameters
        ----------
        forward : mne.Forward
            The mne-python Forward model instance.
        mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
            The MNE data object.
        weight_norm : bool
            Normalize the filter weight matrix W to unit length of the columns.
        alpha : float
            The regularization parameter.

        Return
        ------
        self : object returns itself for convenience

        """
        self.weight_norm = weight_norm
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
        data = self.unpack_data_obj(mne_obj)
        leadfield = self.leadfield
        leadfield /= np.linalg.norm(leadfield, axis=0)
        # leadfield -= leadfield.mean(axis=0)
        n_chans, n_dipoles = self.leadfield.shape

        y = data

        I = np.identity(n_chans)
        y -= y.mean(axis=1, keepdims=True)
        C = self.data_covariance(y, center=False, ddof=1)
        # C = OAS(assume_centered=False).fit(C.T).covariance_.T
        # C = LedoitWolf(assume_centered=False).fit(C.T).covariance_.T

        # Recompute regularization based on the max eigenvalue of the Covariance
        # Matrix (opposed to that of the leadfield)
        # self.alphas = np.logspace(-4, 1, self.n_reg_params) * np.diagonal(y@y.T).mean()
        self.get_alphas(reference=C)

        inverse_operators = []
        for alpha in self.alphas:
            C_inv = self.robust_inverse(C + alpha * I)

            # W = (C_inv @ leadfield) / np.diagonal(leadfield.T @ C_inv @ leadfield)
            upper = C_inv @ leadfield
            lower = np.einsum("ij,jk,ki->i", leadfield.T, C_inv, leadfield)
            W = upper / lower

            # C_inv_L = C_inv @ leadfield
            # diagonal_elements = np.einsum('ij,ji->i', leadfield.T, C_inv_L)
            # W = C_inv_L / diagonal_elements

            if self.weight_norm:
                W /= np.linalg.norm(W, axis=0)

            inverse_operator = W.T
            inverse_operators.append(inverse_operator)

        self.inverse_operators = [
            InverseOperator(inverse_operator, self.name)
            for inverse_operator in inverse_operators
        ]
        return self

__init__ 

__init__(
    name="LCMV Beamformer",
    reduce_rank=True,
    rank="auto",
    **kwargs,
)
Source code in invert/solvers/beamformers/lcmv.py 
37
38
39
def __init__(self, name="LCMV Beamformer", reduce_rank=True, rank="auto", **kwargs):
    self.name = name
    return super().__init__(reduce_rank=reduce_rank, rank=rank, **kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj,
    *args,
    alpha="auto",
    weight_norm=True,
    verbose=0,
    **kwargs,
)

Calculate inverse operator.

Parameters:

Name Type Description Default
forward Forward

The mne-python Forward model instance.

 required
mne_obj [Evoked, Epochs, Raw]

The MNE data object.

 required
weight_norm bool

Normalize the filter weight matrix W to unit length of the columns.

 True
alpha float

The regularization parameter.

 'auto'
Return

self : object returns itself for convenience

Source code in invert/solvers/beamformers/lcmv.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def make_inverse_operator(
    self,
    forward,
    mne_obj,
    *args,
    alpha="auto",
    weight_norm=True,
    verbose=0,
    **kwargs,
):
    """Calculate inverse operator.

    Parameters
    ----------
    forward : mne.Forward
        The mne-python Forward model instance.
    mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
        The MNE data object.
    weight_norm : bool
        Normalize the filter weight matrix W to unit length of the columns.
    alpha : float
        The regularization parameter.

    Return
    ------
    self : object returns itself for convenience

    """
    self.weight_norm = weight_norm
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
    data = self.unpack_data_obj(mne_obj)
    leadfield = self.leadfield
    leadfield /= np.linalg.norm(leadfield, axis=0)
    # leadfield -= leadfield.mean(axis=0)
    n_chans, n_dipoles = self.leadfield.shape

    y = data

    I = np.identity(n_chans)
    y -= y.mean(axis=1, keepdims=True)
    C = self.data_covariance(y, center=False, ddof=1)
    # C = OAS(assume_centered=False).fit(C.T).covariance_.T
    # C = LedoitWolf(assume_centered=False).fit(C.T).covariance_.T

    # Recompute regularization based on the max eigenvalue of the Covariance
    # Matrix (opposed to that of the leadfield)
    # self.alphas = np.logspace(-4, 1, self.n_reg_params) * np.diagonal(y@y.T).mean()
    self.get_alphas(reference=C)

    inverse_operators = []
    for alpha in self.alphas:
        C_inv = self.robust_inverse(C + alpha * I)

        # W = (C_inv @ leadfield) / np.diagonal(leadfield.T @ C_inv @ leadfield)
        upper = C_inv @ leadfield
        lower = np.einsum("ij,jk,ki->i", leadfield.T, C_inv, leadfield)
        W = upper / lower

        # C_inv_L = C_inv @ leadfield
        # diagonal_elements = np.einsum('ij,ji->i', leadfield.T, C_inv_L)
        # W = C_inv_L / diagonal_elements

        if self.weight_norm:
            W /= np.linalg.norm(W, axis=0)

        inverse_operator = W.T
        inverse_operators.append(inverse_operator)

    self.inverse_operators = [
        InverseOperator(inverse_operator, self.name)
        for inverse_operator in inverse_operators
    ]
    return self