Skip to content

ReciPSIICOS (whitened)

Solver ID: RECIPSIICOS_WHITENED

Usage

from invert import Solver

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

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

Overview

Whitened ReciPSIICOS variant that additionally removes correlation subspace components in a whitened representation.

References

  1. Kuznetsova, A., Nurislamova, Y., & Ossadtchi, A. (2021). Modified covariance beamformer for solving MEG inverse problem in the environment with correlated sources. NeuroImage, 228, 117677.

API Reference

Bases: BaseSolver

Whitened ReciPSIICOS Beamformer for M/EEG inverse solution (recommended).

Extends plain ReciPSIICOS by additionally removing correlation subspace components in a whitened representation, providing better suppression of correlated source cross-talk.

References

[1] Kuznetsova, A., Nurislamova, Y., & Ossadtchi, A. (2021). Modified covariance beamformer for solving MEG inverse problem in the environment with correlated sources. Neuroimage, 228, 117677.

Source code in invert/solvers/beamformers/recipsiicos_whitened.py 
 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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
class SolverReciPSIICOSWhitened(BaseSolver):
    """
    Whitened ReciPSIICOS Beamformer for M/EEG inverse solution (recommended).

    Extends plain ReciPSIICOS by additionally removing correlation subspace
    components in a whitened representation, providing better suppression of
    correlated source cross-talk.

    References
    ----------
    [1] Kuznetsova, A., Nurislamova, Y., & Ossadtchi, A. (2021). Modified
        covariance beamformer for solving MEG inverse problem in the environment
        with correlated sources. Neuroimage, 228, 117677.
    """

    meta = SolverMeta(
        slug="recipsiicos_whitened",
        full_name="ReciPSIICOS (whitened)",
        category="Beamformers",
        description=(
            "Whitened ReciPSIICOS variant that additionally removes correlation "
            "subspace components in a whitened representation."
        ),
        references=[
            "Kuznetsova, A., Nurislamova, Y., & Ossadtchi, A. (2021). Modified "
            "covariance beamformer for solving MEG inverse problem in the environment "
            "with correlated sources. NeuroImage, 228, 117677.",
        ],
    )

    def __init__(
        self, name="ReciPSIICOS-Whitened", 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",
        noise_cov=None,
        virtual_sensor_energy=0.99,
        pwr_energy=0.9,
        pwr_rank=None,
        cor_rank=None,
        max_pairs=None,
        seed=0,
        spectral_flip=True,
        verbose=0,
        **kwargs,
    ):
        """
        Calculate inverse operator using whitened ReciPSIICOS method.

        Parameters
        ----------
        forward : mne.Forward
            The mne-python Forward model instance.
        mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
            The MNE data object.
        alpha : float or 'auto'
            The regularization parameter for LCMV backend.
        noise_cov : array, optional
            Noise covariance matrix for whitening.
        virtual_sensor_energy : float
            Fraction of leadfield energy to retain (default: 0.99)
        pwr_energy : float
            Fraction of power subspace energy to retain (default: 0.995)
        pwr_rank : int, optional
            Explicit power subspace rank (overrides pwr_energy)
        cor_rank : int, optional
            Number of correlation modes to remove (default: auto)
        max_pairs : int, optional
            Max source pairs to sample (default: all)
        seed : int
            Random seed for pair sampling (default: 0)
        spectral_flip : bool
            Apply spectral flip to ensure PSD (default: True)

        Returns
        -------
        self : object
            Returns itself for convenience
        """
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
        data = self.unpack_data_obj(mne_obj)

        L = self.leadfield
        Y = data
        n_chans, n_dipoles = L.shape

        # Normalize leadfield columns
        L = L / np.linalg.norm(L, axis=0)

        # 1) Optional noise whitening
        Yw, Lw = _whiten(Y, L, noise_cov)

        # 2) Virtual sensors via leadfield SVD
        UrT, Lr = _virtual_sensors(Lw, keep_energy=virtual_sensor_energy)
        Yr = UrT @ Yw

        # 3) Data covariance in reduced space
        Cr = _cov(Yr)

        # 4) Whitened ReciPSIICOS projection
        Ct = _project_covariance_whitened(
            Cr,
            Lr,
            pwr_energy=pwr_energy,
            pwr_rank=pwr_rank,
            cor_rank=cor_rank,
            max_pairs=max_pairs,
            seed=seed,
        )

        # 5) Spectral flip
        if spectral_flip:
            Ct = _psd_spectral_flip(Ct)

        # 6) Build inverse operators for different regularization values
        # ReciPSIICOS uses trace-shrinkage regularization: C + lam*trace(C)/m*I,
        # where lam is dimensionless. Do NOT eigenvalue-scale it.
        if alpha == "auto":
            reg_values = np.logspace(-6, 1, self.n_reg_params)
        else:
            reg_values = np.asarray([float(alpha)], dtype=float)
        self.alphas = list(np.asarray(reg_values, dtype=float))

        inverse_operators = []
        for reg in self.alphas:
            W_reduced = _lcmv_inverse_operator(Lr, Ct, reg=reg)
            W_full = W_reduced @ UrT
            inverse_operator = W_full
            inverse_operators.append(inverse_operator)

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

__init__ 

__init__(
    name="ReciPSIICOS-Whitened",
    reduce_rank=True,
    rank="auto",
    **kwargs,
)
Source code in invert/solvers/beamformers/recipsiicos_whitened.py 
48
49
50
51
52
def __init__(
    self, name="ReciPSIICOS-Whitened", 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",
    noise_cov=None,
    virtual_sensor_energy=0.99,
    pwr_energy=0.9,
    pwr_rank=None,
    cor_rank=None,
    max_pairs=None,
    seed=0,
    spectral_flip=True,
    verbose=0,
    **kwargs,
)

Calculate inverse operator using whitened ReciPSIICOS method.

Parameters:

Name Type Description Default
forward Forward

The mne-python Forward model instance.

 required
mne_obj [Evoked, Epochs, Raw]

The MNE data object.

 required
alpha float or 'auto'

The regularization parameter for LCMV backend.

 'auto'
noise_cov array

Noise covariance matrix for whitening.

 None
virtual_sensor_energy float

Fraction of leadfield energy to retain (default: 0.99)

 0.99
pwr_energy float

Fraction of power subspace energy to retain (default: 0.995)

 0.9
pwr_rank int

Explicit power subspace rank (overrides pwr_energy)

 None
cor_rank int

Number of correlation modes to remove (default: auto)

 None
max_pairs int

Max source pairs to sample (default: all)

 None
seed int

Random seed for pair sampling (default: 0)

 0
spectral_flip bool

Apply spectral flip to ensure PSD (default: True)

 True

Returns:

Name Type Description
self object

Returns itself for convenience
Source code in invert/solvers/beamformers/recipsiicos_whitened.py 
 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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def make_inverse_operator(
    self,
    forward,
    mne_obj,
    *args,
    alpha="auto",
    noise_cov=None,
    virtual_sensor_energy=0.99,
    pwr_energy=0.9,
    pwr_rank=None,
    cor_rank=None,
    max_pairs=None,
    seed=0,
    spectral_flip=True,
    verbose=0,
    **kwargs,
):
    """
    Calculate inverse operator using whitened ReciPSIICOS method.

    Parameters
    ----------
    forward : mne.Forward
        The mne-python Forward model instance.
    mne_obj : [mne.Evoked, mne.Epochs, mne.io.Raw]
        The MNE data object.
    alpha : float or 'auto'
        The regularization parameter for LCMV backend.
    noise_cov : array, optional
        Noise covariance matrix for whitening.
    virtual_sensor_energy : float
        Fraction of leadfield energy to retain (default: 0.99)
    pwr_energy : float
        Fraction of power subspace energy to retain (default: 0.995)
    pwr_rank : int, optional
        Explicit power subspace rank (overrides pwr_energy)
    cor_rank : int, optional
        Number of correlation modes to remove (default: auto)
    max_pairs : int, optional
        Max source pairs to sample (default: all)
    seed : int
        Random seed for pair sampling (default: 0)
    spectral_flip : bool
        Apply spectral flip to ensure PSD (default: True)

    Returns
    -------
    self : object
        Returns itself for convenience
    """
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
    data = self.unpack_data_obj(mne_obj)

    L = self.leadfield
    Y = data
    n_chans, n_dipoles = L.shape

    # Normalize leadfield columns
    L = L / np.linalg.norm(L, axis=0)

    # 1) Optional noise whitening
    Yw, Lw = _whiten(Y, L, noise_cov)

    # 2) Virtual sensors via leadfield SVD
    UrT, Lr = _virtual_sensors(Lw, keep_energy=virtual_sensor_energy)
    Yr = UrT @ Yw

    # 3) Data covariance in reduced space
    Cr = _cov(Yr)

    # 4) Whitened ReciPSIICOS projection
    Ct = _project_covariance_whitened(
        Cr,
        Lr,
        pwr_energy=pwr_energy,
        pwr_rank=pwr_rank,
        cor_rank=cor_rank,
        max_pairs=max_pairs,
        seed=seed,
    )

    # 5) Spectral flip
    if spectral_flip:
        Ct = _psd_spectral_flip(Ct)

    # 6) Build inverse operators for different regularization values
    # ReciPSIICOS uses trace-shrinkage regularization: C + lam*trace(C)/m*I,
    # where lam is dimensionless. Do NOT eigenvalue-scale it.
    if alpha == "auto":
        reg_values = np.logspace(-6, 1, self.n_reg_params)
    else:
        reg_values = np.asarray([float(alpha)], dtype=float)
    self.alphas = list(np.asarray(reg_values, dtype=float))

    inverse_operators = []
    for reg in self.alphas:
        W_reduced = _lcmv_inverse_operator(Lr, Ct, reg=reg)
        W_full = W_reduced @ UrT
        inverse_operator = W_full
        inverse_operators.append(inverse_operator)

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