ENH: Add Forward-based projection reconstruction - #14235
Conversation
a6d7b36 to
4d386a5
Compare
4d386a5 to
874820f
Compare
larsoner
left a comment
There was a problem hiding this comment.
Looks good to me! Just one idea for supporting more / better rank modes more consistent with other places in the codebase
| if rank is not None: | ||
| inv, _, _ = _reg_pinv(whitened_dots, reg=0, rank=rank) | ||
| elif fmd.get("pinv_method", "tsvd") == "tsvd": | ||
| inv, fmd["nest"] = _pinv_trunc(whitened_dots, fmd["miss"]) |
There was a problem hiding this comment.
This should maybe respect rank actually... maybe we don't need pinv_trunc. (I suspect pinv_trunc is the same as _reg_pinv(..., reg=0, rank=<something>) where the <something> is picked from the data / first singular value.
There was a problem hiding this comment.
yes I had to leave for dinner, did not have time to look at it. will do now.
There was a problem hiding this comment.
Ok I looked into this, in the change I added it was
if rank is not None:
inv, _, _ = _reg_pinv(whitened_dots, reg=0, rank=rank)
the elif fmd.get("pinv_method", "tsvd") == "tsvd": was the original code changed from if to elif. New changes never use _pinv_trunc.
Now for the existing code, I don't think they do the same code and removing _pinv_trunc might require further changes to the function.
|
Just pushed a little simplification, marking for merge-when-green, thanks in advance @BabaSanfour ! |
Reference issue (if any)
Second phase and closes #14208, following #14220.
What does this implement/fix?
Extends
reconstruct_proj()with optionalforwardand integerrankparams.API:
forward=None, rank=None: existing geometry-based reconstruction;forward=Forward, rank=int: Forward-based exact-rank reconstruction;rankwithoutforward, orforwardwithoutrank: for now error. This is what we need in SSP. Acceptingrank/forwardalone would require further tests and more branching.Additional information
Numerical validation
For SSP-SIR (mne-denoise#32, mne-denoise#46), with a projector$P$ and a lead field $L$ , the reconstruction is
which is equivalent in exact arithmetic to$K = LL^\top$ in the sensor-space formulation:
I checked this using TMS-EEG data from three subjects (
CON001,CON006,CON009), 64 EEG channels, and the current MNE-Denoise spherical radial 5000-dipole lead field:1.07e-13at the effective ranks;1.22e-13at the effective ranks;8.37e-14.Limitation
I spent some time playing with where this implementation can fail for a user outside of the intended SSP-SIR goal.
squares the singular spectrum, so Forwards with a very wide singular-value spectrum can lose weak modes in finite precision.
In one strongly ill-conditioned free-orientation EEG Forward used for stress testing:
M=10and became large toward the high-rank tail.AI usage
Ideation and research with GPT-5.6 Sol; implementation assisted by Luna. All tests and code reviewed by me.