Greedy Maximum-Likelihood

Solver ID: GreedyML

Usage

from invert import Solver

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

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

Overview

Source localization via multi-start greedy forward selection with ML objective, BIC model order selection, and coordinate-wise refinement. Scales to large source spaces.

References

1. Wax, M., & Kailath, T. (1985). Detection of signals by information theoretic criteria. IEEE Trans. ASSP, 33(2), 387-392.

API Reference

Bases: BaseSolver

Source localization via multi-start greedy forward selection with ML objective, BIC model order selection, and coordinate-wise refinement.

Scales as O(n_starts * n_sources * k_max) — works for any source space.

References

[1] Wax, M., & Kailath, T. (1985). Detection of signals by information theoretic criteria. IEEE Trans. ASSP, 33(2), 387-392.

Source code in invert/solvers/music/greedy_ml.py

class SolverGreedyML(BaseSolver):
    """Source localization via multi-start greedy forward selection with
    ML objective, BIC model order selection, and coordinate-wise refinement.

    Scales as O(n_starts * n_sources * k_max) — works for any source space.

    References
    ----------
    [1] Wax, M., & Kailath, T. (1985). Detection of signals by information
        theoretic criteria. IEEE Trans. ASSP, 33(2), 387-392.
    """

    meta = SolverMeta(
        acronym="GreedyML",
        full_name="Greedy Maximum-Likelihood",
        category="Subspace Methods",
        description=(
            "Source localization via multi-start greedy forward selection "
            "with ML objective, BIC model order selection, and coordinate-wise "
            "refinement. Scales to large source spaces."
        ),
        references=[
            "Wax, M., & Kailath, T. (1985). Detection of signals by information theoretic criteria. IEEE Trans. ASSP, 33(2), 387-392.",
        ],
    )

    def __init__(self, name="GreedyML", **kwargs):
        self.name = name
        super().__init__(**kwargs)

    def make_inverse_operator(
        self,
        forward,
        mne_obj=None,
        *args,
        alpha="auto",
        noise_cov: mne.Covariance | None = None,
        k_max=5,
        n_starts=50,
        n_refine_iters=3,
        penalty_mode="bic_per_timepoint",
        **kwargs,
    ):
        super().make_inverse_operator(forward, mne_obj, *args, alpha=alpha, **kwargs)
        wf = self.prepare_whitened_forward(noise_cov)

        data = self.unpack_data_obj(mne_obj)
        data_w = wf.sensor_transform @ data
        G_w = wf.G_white

        n_eff, n_sources = G_w.shape
        T = data_w.shape[1]

        Gram = G_w.T @ G_w
        # Ridge for numerical stability (rank-deficient when n_sources >> n_eff)
        eps_gram = 1e-10 * np.trace(Gram) / max(n_sources, 1)
        Gram[np.diag_indices_from(Gram)] += eps_gram
        R = G_w.T @ data_w
        Q = R @ R.T
        total_var = np.sum(data_w**2)

        selected = greedy_ml_search(
            Gram,
            Q,
            total_var,
            n_eff,
            T,
            n_sources,
            k_max,
            n_starts,
            n_refine_iters,
            penalty_mode,
        )
        selected_sources = np.array(selected, dtype=np.intp)

        # Construct ML inverse operator
        G_sel = G_w[:, selected_sources]
        inv_w = np.linalg.lstsq(G_sel, np.eye(n_eff), rcond=None)[0]

        inverse_operator_w = np.zeros((n_sources, n_eff))
        inverse_operator_w[selected_sources, :] = inv_w

        self.inverse_operators = [
            InverseOperator(inverse_operator_w @ wf.sensor_transform, self.name),
        ]
        return self

__init__

__init__(name='GreedyML', **kwargs)

Source code in invert/solvers/music/greedy_ml.py

def __init__(self, name="GreedyML", **kwargs):
    self.name = name
    super().__init__(**kwargs)

make_inverse_operator

make_inverse_operator(
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: Covariance | None = None,
    k_max=5,
    n_starts=50,
    n_refine_iters=3,
    penalty_mode="bic_per_timepoint",
    **kwargs,
)

Source code in invert/solvers/music/greedy_ml.py

def make_inverse_operator(
    self,
    forward,
    mne_obj=None,
    *args,
    alpha="auto",
    noise_cov: mne.Covariance | None = None,
    k_max=5,
    n_starts=50,
    n_refine_iters=3,
    penalty_mode="bic_per_timepoint",
    **kwargs,
):
    super().make_inverse_operator(forward, mne_obj, *args, alpha=alpha, **kwargs)
    wf = self.prepare_whitened_forward(noise_cov)

    data = self.unpack_data_obj(mne_obj)
    data_w = wf.sensor_transform @ data
    G_w = wf.G_white

    n_eff, n_sources = G_w.shape
    T = data_w.shape[1]

    Gram = G_w.T @ G_w
    # Ridge for numerical stability (rank-deficient when n_sources >> n_eff)
    eps_gram = 1e-10 * np.trace(Gram) / max(n_sources, 1)
    Gram[np.diag_indices_from(Gram)] += eps_gram
    R = G_w.T @ data_w
    Q = R @ R.T
    total_var = np.sum(data_w**2)

    selected = greedy_ml_search(
        Gram,
        Q,
        total_var,
        n_eff,
        T,
        n_sources,
        k_max,
        n_starts,
        n_refine_iters,
        penalty_mode,
    )
    selected_sources = np.array(selected, dtype=np.intp)

    # Construct ML inverse operator
    G_sel = G_w[:, selected_sources]
    inv_w = np.linalg.lstsq(G_sel, np.eye(n_eff), rcond=None)[0]

    inverse_operator_w = np.zeros((n_sources, n_eff))
    inverse_operator_w[selected_sources, :] = inv_w

    self.inverse_operators = [
        InverseOperator(inverse_operator_w @ wf.sensor_transform, self.name),
    ]
    return self