Skip to content

Variational Bayes Sparse Bayesian Learning

Solver ID: VB-SBL

Usage

from invert import Solver

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

solver = Solver("VB-SBL")
solver.make_inverse_operator(fwd)
stc = solver.apply_inverse_operator(evoked)
stc.plot()

Overview

Variational Bayesian ARD/SBL solver with Gamma hyperpriors, implemented with efficient sensor-space updates.

References

  1. Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1, 211–244.
  2. Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3), 947–966.

API Reference

Bases: BaseSolver

Variational Bayes Sparse Bayesian Learning (VB-SBL).

Hierarchical model: Y = L X + E X_i,t ~ N(0, γ_i) γ_i governed by a Gamma hyperprior on the corresponding precision.

The implementation uses efficient sensor-space updates (invert m×m matrices) and returns the posterior mean as a linear inverse operator.

Source code in invert/solvers/bayesian/vb_sbl.py 
 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
class SolverVBSBL(BaseSolver):
    """Variational Bayes Sparse Bayesian Learning (VB-SBL).

    Hierarchical model:
        Y = L X + E
        X_i,t ~ N(0, γ_i)
        γ_i governed by a Gamma hyperprior on the corresponding precision.

    The implementation uses efficient sensor-space updates (invert m×m matrices)
    and returns the posterior mean as a linear inverse operator.
    """

    meta = SolverMeta(
        acronym="VB-SBL",
        full_name="Variational Bayes Sparse Bayesian Learning",
        category="Bayesian",
        description=(
            "Variational Bayesian ARD/SBL solver with Gamma hyperpriors, implemented "
            "with efficient sensor-space updates."
        ),
        references=[
            "Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1, 211–244.",
            "Wipf, D., & Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3), 947–966.",
        ],
    )

    def __init__(
        self,
        name: str = "VB-SBL",
        a0: float = 1e-6,
        b0: float = 1e-6,
        **kwargs: Any,
    ) -> None:
        self.name = name
        self.a0 = float(a0)
        self.b0 = float(b0)
        super().__init__(**kwargs)

    def make_inverse_operator(  # type: ignore[override]
        self,
        forward,
        mne_obj,
        *args: Any,
        alpha: str | float = "auto",
        max_iter: int = 300,
        noise_cov: np.ndarray | None = None,
        prune: bool = True,
        pruning_thresh: float = 1e-4,
        convergence_criterion: float = 1e-6,
        **kwargs: Any,
    ):
        super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)

        Y = self.unpack_data_obj(mne_obj)
        Y = Y - Y.mean(axis=0, keepdims=True)

        n_chans = self.leadfield.shape[0]
        if noise_cov is None:
            noise_cov = np.eye(n_chans)

        data_cov = self.data_covariance(Y, center=False, ddof=1)
        self.get_alphas(reference=data_cov)

        inverse_operators = []
        for alpha_eff in self.alphas:
            K = self._vb_sbl(
                Y,
                float(alpha_eff),
                noise_cov=noise_cov,
                max_iter=max_iter,
                prune=prune,
                pruning_thresh=pruning_thresh,
                conv_crit=convergence_criterion,
            )
            inverse_operators.append(K)

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

    def _vb_sbl(
        self,
        Y: np.ndarray,
        noise_var: float,
        *,
        noise_cov: np.ndarray,
        max_iter: int,
        prune: bool,
        pruning_thresh: float,
        conv_crit: float,
    ) -> np.ndarray:
        n_chans, n_times = Y.shape
        L = deepcopy(self.leadfield)
        L_norm = np.linalg.norm(L, axis=0)
        L_norm = np.where(L_norm <= 0, 1.0, L_norm)
        L /= L_norm

        noise_cov = np.asarray(noise_cov, dtype=np.float64)
        noise_cov = noise_cov + 1e-12 * np.eye(noise_cov.shape[0])
        try:
            chol = np.linalg.cholesky(noise_cov)
            Wn = np.linalg.inv(chol)
        except np.linalg.LinAlgError:
            d = np.sqrt(np.clip(np.diag(noise_cov), 1e-12, None))
            Wn = np.diag(1.0 / d)

        Yw = Wn @ Y
        Lw = Wn @ L

        n_dipoles = Lw.shape[1]
        gammas = np.ones(n_dipoles, dtype=np.float64)
        denom_const = self.a0 + 0.5 * float(n_times)

        for _it in range(int(max_iter)):
            Sigma_y = noise_var * np.eye(n_chans) + (Lw * gammas) @ Lw.T
            Sigma_y_inv = self.robust_inverse(Sigma_y)

            L_Sigma = Sigma_y_inv @ Lw
            z_diag = np.sum(Lw * L_Sigma, axis=0)
            diag_Sigma_x = gammas - gammas * gammas * z_diag

            mu_x = (gammas[:, None] * (Lw.T @ (Sigma_y_inv @ Yw))).astype(np.float64)

            Ex2 = np.sum(mu_x * mu_x, axis=1) + float(n_times) * diag_Sigma_x
            gammas_new = (self.b0 + 0.5 * Ex2) / denom_const

            if prune:
                thr = float(pruning_thresh) * float(np.max(gammas_new))
                gammas_new = np.where(gammas_new >= thr, gammas_new, 0.0)

            rel = float(np.linalg.norm(gammas_new - gammas)) / max(
                float(np.linalg.norm(gammas)), 1e-15
            )
            gammas = gammas_new
            if rel < float(conv_crit):
                break

        Sigma_y = noise_var * np.eye(n_chans) + (Lw * gammas) @ Lw.T
        Sigma_y_inv = self.robust_inverse(Sigma_y)

        temp = Sigma_y_inv @ Wn
        K_scaled = gammas[:, None] * (Lw.T @ temp)
        K = (K_scaled / L_norm[:, None]).astype(np.float64)
        return K

__init__ 

__init__(
    name: str = "VB-SBL",
    a0: float = 1e-06,
    b0: float = 1e-06,
    **kwargs: Any,
) -> None
Source code in invert/solvers/bayesian/vb_sbl.py 
40
41
42
43
44
45
46
47
48
49
50
def __init__(
    self,
    name: str = "VB-SBL",
    a0: float = 1e-6,
    b0: float = 1e-6,
    **kwargs: Any,
) -> None:
    self.name = name
    self.a0 = float(a0)
    self.b0 = float(b0)
    super().__init__(**kwargs)

make_inverse_operator 

make_inverse_operator(
    forward,
    mne_obj,
    *args: Any,
    alpha: str | float = "auto",
    max_iter: int = 300,
    noise_cov: ndarray | None = None,
    prune: bool = True,
    pruning_thresh: float = 0.0001,
    convergence_criterion: float = 1e-06,
    **kwargs: Any,
)
Source code in invert/solvers/bayesian/vb_sbl.py 
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
def make_inverse_operator(  # type: ignore[override]
    self,
    forward,
    mne_obj,
    *args: Any,
    alpha: str | float = "auto",
    max_iter: int = 300,
    noise_cov: np.ndarray | None = None,
    prune: bool = True,
    pruning_thresh: float = 1e-4,
    convergence_criterion: float = 1e-6,
    **kwargs: Any,
):
    super().make_inverse_operator(forward, *args, alpha=alpha, **kwargs)

    Y = self.unpack_data_obj(mne_obj)
    Y = Y - Y.mean(axis=0, keepdims=True)

    n_chans = self.leadfield.shape[0]
    if noise_cov is None:
        noise_cov = np.eye(n_chans)

    data_cov = self.data_covariance(Y, center=False, ddof=1)
    self.get_alphas(reference=data_cov)

    inverse_operators = []
    for alpha_eff in self.alphas:
        K = self._vb_sbl(
            Y,
            float(alpha_eff),
            noise_cov=noise_cov,
            max_iter=max_iter,
            prune=prune,
            pruning_thresh=pruning_thresh,
            conv_crit=convergence_criterion,
        )
        inverse_operators.append(K)

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