Skip to content

Empirical Bayesian Beamformer

Solver ID: EBB

Usage

from invert import Solver

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

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

Overview

Iterative empirical-Bayes / ARD-style beamformer that updates per-source variances (hyperparameters) from the data and forms corresponding spatial filters.

References

  1. tbd

API Reference

Bases: BaseSolver

Empirical Bayesian Beamformer (EBB) solver for M/EEG inverse problem.

Source code in invert/solvers/beamformers/ebb.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
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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
class SolverEBB(BaseSolver):
    """
    Empirical Bayesian Beamformer (EBB) solver for M/EEG inverse problem.
    """

    meta = SolverMeta(
        slug="ebb",
        full_name="Empirical Bayesian Beamformer",
        category="Beamformers",
        description=(
            "Iterative empirical-Bayes / ARD-style beamformer that updates per-source "
            "variances (hyperparameters) from the data and forms corresponding spatial "
            "filters."
        ),
        references=["tbd"],
    )

    def __init__(
        self,
        name="Empirical Bayesian 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,
        weight_norm=True,
        noise_cov=None,
        alpha="auto",
        max_iter=100,
        tol=1e-6,
        **kwargs,
    ):
        """
        Solve the inverse problem using the Empirical Bayesian Beamformer method.

        This implements an iterative algorithm that estimates source variances (hyperparameters)
        from the data using an empirical Bayes approach, similar to Champagne/ARD methods.

        Parameters:
        -----------
        data : array, shape (n_channels, n_times)
            The sensor data.
        forward : array, shape (n_channels, n_sources)
            The forward solution.
        noise_cov : array, shape (n_channels, n_channels), optional
            The noise covariance matrix.
        max_iter : int
            Maximum number of iterations for the EM algorithm.
        tol : float
            Convergence tolerance for source variance updates.

        Returns:
        --------
        sources : array, shape (n_sources, n_times)
            The estimated source time series.
        """
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)

        data = self.unpack_data_obj(mne_obj)
        leadfield = self.leadfield
        n_channels, n_times = data.shape
        n_sources = leadfield.shape[1]

        # normalize leadfield
        leadfield /= np.linalg.norm(leadfield, axis=0)

        # Compute data covariance
        data_cov = self.data_covariance(data, center=True, ddof=1)

        # handle noise_cov - use scaled identity as default
        if noise_cov is None:
            noise_cov = np.identity(n_channels)

        inverse_operators = []
        self.alphas = self.get_alphas(reference=data_cov)

        for alpha in self.alphas:
            # Initialize source variances (diagonal of source covariance)
            # Start with uniform prior
            gamma = np.ones(n_sources)

            # Regularized noise covariance
            Cn = alpha * noise_cov

            # Iterative empirical Bayes updates
            for n_iter in range(max_iter):
                logger.debug(f"EBB iteration {n_iter + 1} of {max_iter}")
                # Build current model covariance: C = L @ diag(gamma) @ L.T + Cn
                # For efficiency, compute using woodbury: C^-1 = Cn^-1 - Cn^-1 @ L @ (I/gamma + L.T @ Cn^-1 @ L)^-1 @ L.T @ Cn^-1

                Cn_inv = self.robust_inverse(Cn)

                # Compute L.T @ Cn^-1 @ L efficiently
                Cn_inv_L = Cn_inv @ leadfield
                LT_Cn_inv_L = leadfield.T @ Cn_inv_L

                # Add diagonal loading: (I/gamma + L.T @ Cn^-1 @ L)
                middle = np.diag(1.0 / gamma) + LT_Cn_inv_L
                middle_inv = self.robust_inverse(middle)

                # Woodbury identity for C^-1
                Cn_inv - Cn_inv_L @ middle_inv @ Cn_inv_L.T

                # Compute posterior source covariance: Sigma_s = (I/gamma + L.T @ C^-1 @ L)^-1
                # This is actually just middle_inv from above!
                Sigma_s = middle_inv

                # Compute beamformer weights for this iteration
                W = Sigma_s @ leadfield.T @ Cn_inv  # shape: (n_sources, n_channels)

                # Estimate source activity
                source_estimates = W @ data  # shape: (n_sources, n_times)

                # Update source variances using empirical estimates
                # gamma_new[i] = trace(Sigma_s[i,i]) + (source_estimates[i,:]^2).mean()
                # Simplified: use source power + trace correction
                source_power = np.mean(
                    source_estimates**2, axis=1
                )  # shape: (n_sources,)
                gamma_new = source_power + np.diag(Sigma_s)

                # Check convergence
                rel_change = np.abs(gamma_new - gamma) / (np.abs(gamma) + 1e-10)
                max_change = np.max(rel_change)

                if max_change < tol:
                    logger.info(
                        f"EBB converged after {n_iter + 1} iterations (max change: {max_change:.2e})"
                    )
                    break

                # Update with damping for stability
                damping = 0.5
                gamma = damping * gamma_new + (1 - damping) * gamma

                # Prevent numerical issues
                gamma = np.maximum(gamma, 1e-12)

            else:
                logger.warning(
                    f"EBB reached max iterations ({max_iter}), max change: {max_change:.2e}"
                )

            # Final beamformer weights with converged hyperparameters
            Cn_inv = self.robust_inverse(Cn)
            Cn_inv_L = Cn_inv @ leadfield
            LT_Cn_inv_L = leadfield.T @ Cn_inv_L
            middle = np.diag(1.0 / gamma) + LT_Cn_inv_L
            middle_inv = self.robust_inverse(middle)
            Cn_inv - Cn_inv_L @ middle_inv @ Cn_inv_L.T

            W = middle_inv @ leadfield.T @ Cn_inv

            # Optional weight normalization
            if weight_norm:
                W /= np.linalg.norm(W, axis=1, keepdims=True)

            inverse_operators.append(W)  # Transpose to match expected shape

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

        return self

__init__ 

__init__(
    name="Empirical Bayesian Beamformer",
    reduce_rank=True,
    rank="auto",
    **kwargs,
)
Source code in invert/solvers/beamformers/ebb.py 
27
28
29
30
31
32
33
34
35
def __init__(
    self,
    name="Empirical Bayesian 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,
    weight_norm=True,
    noise_cov=None,
    alpha="auto",
    max_iter=100,
    tol=1e-06,
    **kwargs,
)

Solve the inverse problem using the Empirical Bayesian Beamformer method.

This implements an iterative algorithm that estimates source variances (hyperparameters) from the data using an empirical Bayes approach, similar to Champagne/ARD methods.

Parameters:

data : array, shape (n_channels, n_times) The sensor data. forward : array, shape (n_channels, n_sources) The forward solution. noise_cov : array, shape (n_channels, n_channels), optional The noise covariance matrix. max_iter : int Maximum number of iterations for the EM algorithm. tol : float Convergence tolerance for source variance updates.

Returns:

sources : array, shape (n_sources, n_times) The estimated source time series.

Source code in invert/solvers/beamformers/ebb.py 
 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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def make_inverse_operator(
    self,
    forward,
    mne_obj,
    *args,
    weight_norm=True,
    noise_cov=None,
    alpha="auto",
    max_iter=100,
    tol=1e-6,
    **kwargs,
):
    """
    Solve the inverse problem using the Empirical Bayesian Beamformer method.

    This implements an iterative algorithm that estimates source variances (hyperparameters)
    from the data using an empirical Bayes approach, similar to Champagne/ARD methods.

    Parameters:
    -----------
    data : array, shape (n_channels, n_times)
        The sensor data.
    forward : array, shape (n_channels, n_sources)
        The forward solution.
    noise_cov : array, shape (n_channels, n_channels), optional
        The noise covariance matrix.
    max_iter : int
        Maximum number of iterations for the EM algorithm.
    tol : float
        Convergence tolerance for source variance updates.

    Returns:
    --------
    sources : array, shape (n_sources, n_times)
        The estimated source time series.
    """
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)

    data = self.unpack_data_obj(mne_obj)
    leadfield = self.leadfield
    n_channels, n_times = data.shape
    n_sources = leadfield.shape[1]

    # normalize leadfield
    leadfield /= np.linalg.norm(leadfield, axis=0)

    # Compute data covariance
    data_cov = self.data_covariance(data, center=True, ddof=1)

    # handle noise_cov - use scaled identity as default
    if noise_cov is None:
        noise_cov = np.identity(n_channels)

    inverse_operators = []
    self.alphas = self.get_alphas(reference=data_cov)

    for alpha in self.alphas:
        # Initialize source variances (diagonal of source covariance)
        # Start with uniform prior
        gamma = np.ones(n_sources)

        # Regularized noise covariance
        Cn = alpha * noise_cov

        # Iterative empirical Bayes updates
        for n_iter in range(max_iter):
            logger.debug(f"EBB iteration {n_iter + 1} of {max_iter}")
            # Build current model covariance: C = L @ diag(gamma) @ L.T + Cn
            # For efficiency, compute using woodbury: C^-1 = Cn^-1 - Cn^-1 @ L @ (I/gamma + L.T @ Cn^-1 @ L)^-1 @ L.T @ Cn^-1

            Cn_inv = self.robust_inverse(Cn)

            # Compute L.T @ Cn^-1 @ L efficiently
            Cn_inv_L = Cn_inv @ leadfield
            LT_Cn_inv_L = leadfield.T @ Cn_inv_L

            # Add diagonal loading: (I/gamma + L.T @ Cn^-1 @ L)
            middle = np.diag(1.0 / gamma) + LT_Cn_inv_L
            middle_inv = self.robust_inverse(middle)

            # Woodbury identity for C^-1
            Cn_inv - Cn_inv_L @ middle_inv @ Cn_inv_L.T

            # Compute posterior source covariance: Sigma_s = (I/gamma + L.T @ C^-1 @ L)^-1
            # This is actually just middle_inv from above!
            Sigma_s = middle_inv

            # Compute beamformer weights for this iteration
            W = Sigma_s @ leadfield.T @ Cn_inv  # shape: (n_sources, n_channels)

            # Estimate source activity
            source_estimates = W @ data  # shape: (n_sources, n_times)

            # Update source variances using empirical estimates
            # gamma_new[i] = trace(Sigma_s[i,i]) + (source_estimates[i,:]^2).mean()
            # Simplified: use source power + trace correction
            source_power = np.mean(
                source_estimates**2, axis=1
            )  # shape: (n_sources,)
            gamma_new = source_power + np.diag(Sigma_s)

            # Check convergence
            rel_change = np.abs(gamma_new - gamma) / (np.abs(gamma) + 1e-10)
            max_change = np.max(rel_change)

            if max_change < tol:
                logger.info(
                    f"EBB converged after {n_iter + 1} iterations (max change: {max_change:.2e})"
                )
                break

            # Update with damping for stability
            damping = 0.5
            gamma = damping * gamma_new + (1 - damping) * gamma

            # Prevent numerical issues
            gamma = np.maximum(gamma, 1e-12)

        else:
            logger.warning(
                f"EBB reached max iterations ({max_iter}), max change: {max_change:.2e}"
            )

        # Final beamformer weights with converged hyperparameters
        Cn_inv = self.robust_inverse(Cn)
        Cn_inv_L = Cn_inv @ leadfield
        LT_Cn_inv_L = leadfield.T @ Cn_inv_L
        middle = np.diag(1.0 / gamma) + LT_Cn_inv_L
        middle_inv = self.robust_inverse(middle)
        Cn_inv - Cn_inv_L @ middle_inv @ Cn_inv_L.T

        W = middle_inv @ leadfield.T @ Cn_inv

        # Optional weight normalization
        if weight_norm:
            W /= np.linalg.norm(W, axis=1, keepdims=True)

        inverse_operators.append(W)  # Transpose to match expected shape

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

    return self