Mixed-Norm Estimate (BCD with Active Set)¶

Solver ID: MxNE-BCD

Usage¶

from invert import Solver

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

solver = Solver("MxNE-BCD")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview¶

L21 mixed-norm inverse via block coordinate descent with active-set screening and optional amplitude debiasing.

References¶

Gramfort, A., Kowalski, M., & Hämäläinen, M. (2012). Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods. Physics in Medicine and Biology, 57(7), 1937-1961.

API Reference¶

Bases: BaseSolver

Mixed-Norm Estimate via Block Coordinate Descent with Active Set.

Solves the L21 mixed-norm problem::

min_X  0.5 ||M - G X||_F^2  +  alpha * sum_j ||X_j||_F

The L21 penalty enforces spatial sparsity (few active source positions) while each active source retains a free time course. The BCD solver with active-set screening avoids processing the full source space, giving large speedups when the true source configuration is sparse.

References

Gramfort, A., Kowalski, M., & Hämäläinen, M. (2012). Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods. Physics in Medicine and Biology, 57(7), 1937-1961.

Source code in invert/solvers/minimum_norm/mxne.py

class SolverMxNEBCD(BaseSolver):
    """Mixed-Norm Estimate via Block Coordinate Descent with Active Set.

    Solves the L21 mixed-norm problem::

        min_X  0.5 ||M - G X||_F^2  +  alpha * sum_j ||X_j||_F

    The L21 penalty enforces spatial sparsity (few active source positions)
    while each active source retains a free time course.  The BCD solver
    with active-set screening avoids processing the full source space,
    giving large speedups when the true source configuration is sparse.

    References
    ----------
    Gramfort, A., Kowalski, M., & Hämäläinen, M. (2012). Mixed-norm
    estimates for the M/EEG inverse problem using accelerated gradient
    methods. Physics in Medicine and Biology, 57(7), 1937-1961.
    """

    meta = SolverMeta(
        acronym="MxNE-BCD",
        full_name="Mixed-Norm Estimate (BCD with Active Set)",
        category="Sparse / Mixed Norm",
        description=(
            "L21 mixed-norm inverse via block coordinate descent with "
            "active-set screening and optional amplitude debiasing."
        ),
        references=[
            "Gramfort, A., Kowalski, M., & Hämäläinen, M. (2012). Mixed-norm "
            "estimates for the M/EEG inverse problem using accelerated gradient "
            "methods. Physics in Medicine and Biology, 57(7), 1937-1961.",
        ],
    )

    SUPPORTS_VECTOR_ORIENTATION = True

    def __init__(self, name="MxNE-BCD", **kwargs):
        self.name = name
        super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        *args,
        alpha=0.01,
        noise_cov: mne.Covariance | None = None,
        **kwargs,
    ):
        """Prepare whitened forward model.

        Parameters
        ----------
        forward : mne.Forward
        alpha : float
            Unused here (regularization is set in apply_inverse_operator).
        noise_cov : mne.Covariance | None
            Noise covariance for sensor whitening.
        """
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
        self.prepare_whitened_forward(noise_cov)
        return self

    def apply_inverse_operator(
        self,
        mne_obj,
        alpha=0.2,
        maxit=3000,
        tol=1e-4,
        active_set_size=10,
        dgap_freq=10,
        bcd_maxit=200,
        debias=True,
        time_pca=True,
        center_data=True,
        depth=0.8,
        depth_limit=10.0,
    ) -> mne.SourceEstimate:
        """Apply the MxNE-BCD solver.

        Parameters
        ----------
        mne_obj : mne.Evoked | mne.Epochs | mne.io.Raw
        alpha : float
            Regularization as a fraction of alpha_max (0 = no penalty,
            1 = all-zero solution).  Typical range 0.05 – 0.5.
        maxit : int
            Maximum outer (active-set expansion) iterations.
        tol : float
            Relative duality-gap tolerance for convergence.
        active_set_size : int
            Number of source positions added per expansion step.
        dgap_freq : int
            Evaluate duality gap every *dgap_freq* BCD iterations.
        bcd_maxit : int
            Maximum BCD iterations per active-set solve.
        debias : bool
            Apply FISTA debiasing to correct L1 amplitude shrinkage.
        time_pca : bool
            Reduce temporal dimension via SVD when n_times > n_sensors.
        center_data : bool
            Subtract per-channel temporal mean before solving.
        depth : float
            Depth weighting exponent applied to leadfield column norms.
            Equalizes sensitivity across source depths (0 = off, 1 = full
            column normalization, 0.8 = MNE default).
        depth_limit : float
            Maximum ratio between largest and smallest depth weight.
        """
        data = self.unpack_data_obj(mne_obj)
        self.validate_operator_data_compatibility(data)
        data = self._sensor_transform @ data

        G = self.leadfield.copy()
        n_orient = int(getattr(self, "_n_orient", 1))
        n_sensors, n_sources = G.shape
        n_positions = n_sources // n_orient

        M = data.copy()
        if center_data:
            M -= M.mean(axis=0, keepdims=True)

        # ------------------------------------------------------------------
        # Depth weighting: normalize leadfield columns to equalize
        # sensitivity across source depths.
        # ------------------------------------------------------------------
        depth_w = None
        if depth > 0:
            depth_w = _compute_depth_weights(
                G,
                n_orient,
                depth,
                depth_limit,
            )
            G = G / depth_w[np.newaxis, :]

        # Time PCA: compress temporal dimension for efficiency
        Vh = None
        if time_pca and M.shape[1] > n_sensors:
            U, s, Vh = np.linalg.svd(M, full_matrices=False)
            rank = max(1, int(np.sum(s > s[0] * 1e-6)))
            M = U[:, :rank] * s[:rank]
            Vh = Vh[:rank]

        n_times = M.shape[1]
        if n_times == 0 or np.allclose(M, 0):
            return self.source_to_object(np.zeros((n_sources, data.shape[1])))

        # ------------------------------------------------------------------
        # Regularization: alpha as fraction of alpha_max
        # ------------------------------------------------------------------
        GtM = G.T @ M
        alpha_max = _norm_l2inf(GtM, n_orient)
        if alpha_max == 0:
            return self.source_to_object(np.zeros((n_sources, data.shape[1])))
        alpha_eff = float(alpha) * alpha_max

        # ------------------------------------------------------------------
        # Lipschitz constants (precomputed, reused across iterations)
        # ------------------------------------------------------------------
        lipschitz = _compute_lipschitz(G, n_orient)

        # ------------------------------------------------------------------
        # Initialise solution and active set
        # ------------------------------------------------------------------
        X = np.zeros((n_sources, n_times))
        R = M.copy()

        gn2 = _groups_norm2(GtM, n_orient)
        n_init = min(active_set_size, n_positions)
        active = set(np.argsort(gn2)[-n_init:].tolist())

        # ------------------------------------------------------------------
        # Outer loop: BCD solve → expand active set → repeat
        # ------------------------------------------------------------------
        converged = False
        for _outer in range(maxit):
            # --- BCD iterations on current active set ---
            for bcd_it in range(bcd_maxit):
                active = _bcd_pass(
                    G,
                    X,
                    R,
                    sorted(active),
                    lipschitz,
                    n_orient,
                    alpha_eff,
                )

                if (bcd_it + 1) % dgap_freq == 0:
                    gap, p_obj = _dgap_l21(G, X, R, alpha_eff, n_orient)
                    if gap < tol * max(p_obj, 1e-30):
                        converged = True
                        break

            if converged:
                break

            # --- Expand active set ---
            GtR = G.T @ R
            gn2_resid = _groups_norm2(GtR, n_orient)
            for j in active:
                gn2_resid[j] = 0.0

            remaining = n_positions - len(active)
            n_add = min(active_set_size, remaining)
            if n_add <= 0:
                break

            new_pos = np.argsort(gn2_resid)[-n_add:]
            active |= set(new_pos.tolist())

        # ------------------------------------------------------------------
        # Optional debiasing
        # ------------------------------------------------------------------
        if debias:
            X = _debias(G, M, X, n_orient)

        # ------------------------------------------------------------------
        # Undo depth weighting and time PCA
        # ------------------------------------------------------------------
        if depth_w is not None:
            X = X / depth_w[:, np.newaxis]

        if Vh is not None:
            X = X @ Vh

        return self.source_to_object(X)

init ¶

__init__(name='MxNE-BCD', **kwargs)

Source code in invert/solvers/minimum_norm/mxne.py

def __init__(self, name="MxNE-BCD", **kwargs):
    self.name = name
    super().__init__(**kwargs)

make_inverse_operator ¶

make_inverse_operator(
    forward,
    *args,
    alpha=0.01,
    noise_cov: Covariance | None = None,
    **kwargs,
)

Prepare whitened forward model.

Parameters:

Name	Type	Description	Default
`forward`	`Forward`		required
`alpha`	`float`	Unused here (regularization is set in apply_inverse_operator).	`0.01`
`noise_cov`	`Covariance \| None`	Noise covariance for sensor whitening.	`None`

Source code in invert/solvers/minimum_norm/mxne.py

def make_inverse_operator(
    self,
    forward,
    *args,
    alpha=0.01,
    noise_cov: mne.Covariance | None = None,
    **kwargs,
):
    """Prepare whitened forward model.

    Parameters
    ----------
    forward : mne.Forward
    alpha : float
        Unused here (regularization is set in apply_inverse_operator).
    noise_cov : mne.Covariance | None
        Noise covariance for sensor whitening.
    """
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)
    self.prepare_whitened_forward(noise_cov)
    return self

apply_inverse_operator ¶

apply_inverse_operator(
    mne_obj,
    alpha=0.2,
    maxit=3000,
    tol=0.0001,
    active_set_size=10,
    dgap_freq=10,
    bcd_maxit=200,
    debias=True,
    time_pca=True,
    center_data=True,
    depth=0.8,
    depth_limit=10.0,
) -> mne.SourceEstimate

Apply the MxNE-BCD solver.

Parameters:

Name	Type	Description	Default
`mne_obj`	`Evoked \| Epochs \| Raw`		required
`alpha`	`float`	Regularization as a fraction of alpha_max (0 = no penalty, 1 = all-zero solution). Typical range 0.05 – 0.5.	`0.2`
`maxit`	`int`	Maximum outer (active-set expansion) iterations.	`3000`
`tol`	`float`	Relative duality-gap tolerance for convergence.	`0.0001`
`active_set_size`	`int`	Number of source positions added per expansion step.	`10`
`dgap_freq`	`int`	Evaluate duality gap every dgap_freq BCD iterations.	`10`
`bcd_maxit`	`int`	Maximum BCD iterations per active-set solve.	`200`
`debias`	`bool`	Apply FISTA debiasing to correct L1 amplitude shrinkage.	`True`
`time_pca`	`bool`	Reduce temporal dimension via SVD when n_times > n_sensors.	`True`
`center_data`	`bool`	Subtract per-channel temporal mean before solving.	`True`
`depth`	`float`	Depth weighting exponent applied to leadfield column norms. Equalizes sensitivity across source depths (0 = off, 1 = full column normalization, 0.8 = MNE default).	`0.8`
`depth_limit`	`float`	Maximum ratio between largest and smallest depth weight.	`10.0`

Source code in invert/solvers/minimum_norm/mxne.py

def apply_inverse_operator(
    self,
    mne_obj,
    alpha=0.2,
    maxit=3000,
    tol=1e-4,
    active_set_size=10,
    dgap_freq=10,
    bcd_maxit=200,
    debias=True,
    time_pca=True,
    center_data=True,
    depth=0.8,
    depth_limit=10.0,
) -> mne.SourceEstimate:
    """Apply the MxNE-BCD solver.

    Parameters
    ----------
    mne_obj : mne.Evoked | mne.Epochs | mne.io.Raw
    alpha : float
        Regularization as a fraction of alpha_max (0 = no penalty,
        1 = all-zero solution).  Typical range 0.05 – 0.5.
    maxit : int
        Maximum outer (active-set expansion) iterations.
    tol : float
        Relative duality-gap tolerance for convergence.
    active_set_size : int
        Number of source positions added per expansion step.
    dgap_freq : int
        Evaluate duality gap every *dgap_freq* BCD iterations.
    bcd_maxit : int
        Maximum BCD iterations per active-set solve.
    debias : bool
        Apply FISTA debiasing to correct L1 amplitude shrinkage.
    time_pca : bool
        Reduce temporal dimension via SVD when n_times > n_sensors.
    center_data : bool
        Subtract per-channel temporal mean before solving.
    depth : float
        Depth weighting exponent applied to leadfield column norms.
        Equalizes sensitivity across source depths (0 = off, 1 = full
        column normalization, 0.8 = MNE default).
    depth_limit : float
        Maximum ratio between largest and smallest depth weight.
    """
    data = self.unpack_data_obj(mne_obj)
    self.validate_operator_data_compatibility(data)
    data = self._sensor_transform @ data

    G = self.leadfield.copy()
    n_orient = int(getattr(self, "_n_orient", 1))
    n_sensors, n_sources = G.shape
    n_positions = n_sources // n_orient

    M = data.copy()
    if center_data:
        M -= M.mean(axis=0, keepdims=True)

    # ------------------------------------------------------------------
    # Depth weighting: normalize leadfield columns to equalize
    # sensitivity across source depths.
    # ------------------------------------------------------------------
    depth_w = None
    if depth > 0:
        depth_w = _compute_depth_weights(
            G,
            n_orient,
            depth,
            depth_limit,
        )
        G = G / depth_w[np.newaxis, :]

    # Time PCA: compress temporal dimension for efficiency
    Vh = None
    if time_pca and M.shape[1] > n_sensors:
        U, s, Vh = np.linalg.svd(M, full_matrices=False)
        rank = max(1, int(np.sum(s > s[0] * 1e-6)))
        M = U[:, :rank] * s[:rank]
        Vh = Vh[:rank]

    n_times = M.shape[1]
    if n_times == 0 or np.allclose(M, 0):
        return self.source_to_object(np.zeros((n_sources, data.shape[1])))

    # ------------------------------------------------------------------
    # Regularization: alpha as fraction of alpha_max
    # ------------------------------------------------------------------
    GtM = G.T @ M
    alpha_max = _norm_l2inf(GtM, n_orient)
    if alpha_max == 0:
        return self.source_to_object(np.zeros((n_sources, data.shape[1])))
    alpha_eff = float(alpha) * alpha_max

    # ------------------------------------------------------------------
    # Lipschitz constants (precomputed, reused across iterations)
    # ------------------------------------------------------------------
    lipschitz = _compute_lipschitz(G, n_orient)

    # ------------------------------------------------------------------
    # Initialise solution and active set
    # ------------------------------------------------------------------
    X = np.zeros((n_sources, n_times))
    R = M.copy()

    gn2 = _groups_norm2(GtM, n_orient)
    n_init = min(active_set_size, n_positions)
    active = set(np.argsort(gn2)[-n_init:].tolist())

    # ------------------------------------------------------------------
    # Outer loop: BCD solve → expand active set → repeat
    # ------------------------------------------------------------------
    converged = False
    for _outer in range(maxit):
        # --- BCD iterations on current active set ---
        for bcd_it in range(bcd_maxit):
            active = _bcd_pass(
                G,
                X,
                R,
                sorted(active),
                lipschitz,
                n_orient,
                alpha_eff,
            )

            if (bcd_it + 1) % dgap_freq == 0:
                gap, p_obj = _dgap_l21(G, X, R, alpha_eff, n_orient)
                if gap < tol * max(p_obj, 1e-30):
                    converged = True
                    break

        if converged:
            break

        # --- Expand active set ---
        GtR = G.T @ R
        gn2_resid = _groups_norm2(GtR, n_orient)
        for j in active:
            gn2_resid[j] = 0.0

        remaining = n_positions - len(active)
        n_add = min(active_set_size, remaining)
        if n_add <= 0:
            break

        new_pos = np.argsort(gn2_resid)[-n_add:]
        active |= set(new_pos.tolist())

    # ------------------------------------------------------------------
    # Optional debiasing
    # ------------------------------------------------------------------
    if debias:
        X = _debias(G, M, X, n_orient)

    # ------------------------------------------------------------------
    # Undo depth weighting and time PCA
    # ------------------------------------------------------------------
    if depth_w is not None:
        X = X / depth_w[:, np.newaxis]

    if Vh is not None:
        X = X @ Vh

    return self.source_to_object(X)

Mixed-Norm Estimate (BCD with Active Set)¶

Usage¶

Overview¶

References¶

API Reference¶

__init__ ¶

make_inverse_operator ¶

apply_inverse_operator ¶

init ¶