Add noise-weighted Maxwell Filtering (SSS) using Foster's Inverse

Start by reading mne.compute_raw_covariance, mne.preprocessing.compute_maxwell_basis, and mne.preprocessing.maxwell_filter to trace the existing SSS calculation and multipole-moment pseudoinverse. The work is done when Maxwell Filtering supports the proposed noise-weighted Foster's Inverse using the covariance options described here, with results demonstrating more stable, lower-noise OPM-MEG reconstruction.

Describe the new feature or enhancement

MEG data is often preprocessed with Maxwell Filtering and the Signal Space Separation method (SSS), which separates the brain activity from other measured magnetic fields. Under conditions of high sensor noise, the SSS method can become unstable. We propose a novel application of a weighted matrix inverse method for the estimate of each magnetic topography order (multipole moment weights), called Foster's Inverse, where the SSS calculation is modified to take the sensor noise covariance into account. This method is developed by myself, @larsoner, and Dr. Samu Taulu and helps mitigate higher sensor noise and sensitivity to low-frequency magnetic fields seen with OPM-MEG, specifically. This method is in the process of publication, and we would like to add the feature to MNE-Python.

Describe your proposed implementation

The noise covariance matrix N is calculated using the Empirical method already implemented in MNE-Python from the raw data baseline period using:
N = mne.compute_raw_covariance(raw,rank="info",method='empirical')["data"]

From here, the SSS basis calculations proceed as usual with mne.preprocessing.compute_maxwell_basis, which gives the spatial matrix S and multipole moments reg_moments. The typical next step within mne.preprocessing.maxwell_filter is to estimate the multipole moments by multiplying the data matrix with the psuedoinverse of S, pS.

This step can be unstable with high sensor noise, causing the spreading of artifacts across channels. We propose to modify this step with Foster's Inverse, where the main extra step is modifying the pseudoinverse to be weighted with matrix N:
np.linalg.pinv(S@A@S_star +N)
where A is the normalized covariance of the initial estimate for the multipole moments reg_moments.

Describe possible alternatives

N can be calculated in other ways, such as from an empty-room recording, or using methods other than Empirical. Our work suggests that the empirical method is the simplest and functions well. Also, the matrix S can be calculated using other means, such as the multi-origin method proposed by McPherson et al., or spheroidal harmonics proposed by Tierney et al. Both could be explored in the future. Our proposed method is a simple and robust adaptation of already existing methodology within MNE-Python. Preliminary results with OPM-MEG data show stable reconstructions of brain activity with reduced noise and improved source localization.

Additional context

We don't have a preprint yet, but the Foster's Inverse with SSS method and preliminary results are discussed in Chapter 4 of my PhD dissertation which can be found here: Improving the methodology and instrumentation for on-scalp magnetoencephalography (MEG)