Size of resolution matrix (make_inverse_resolution_matrix) inconsistent with documentation

[Katharinski]

Describe the bug

It is unclear what the correct way to compute the resolution matrix is when using mne.minimum_norm.make_inverse_resolution_matrix with an inverse operator that has loose orientations. The function _get_matrix_from_inverse_operator is called (in line 59) to return a matrix that can be multiplied with the leadfield, and the documentation there states:

For loose orientation constraint, the CTFs are computed for the normal component (pick_ori=‘normal’).

However, it does not seem like pick_ori='normal' is actually used.

Steps to reproduce

import mne 

data_path = mne.datasets.sample.data_path()
fname_inv = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-eeg-inv.fif' # free source orientation
inverse_operator = mne.minimum_norm.read_inverse_operator(fname_inv)

fname_fwd = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif'
forward = mne.read_forward_solution(fname_fwd)

method = "sLORETA" 
snr = 3.
lambda2 = 1. / snr ** 2

res_matrix = mne.minimum_norm.make_inverse_resolution_matrix(forward,inverse_operator,method=method,lambda2=lambda2)

Expected results

I expect the res_matrix to have dimensions nsources x nsources, where nsources is the number of sources in the source space.

Actual results

The res_matrix has dimensions (3nsources) x (3nsources).

Additional information

Output of mne.sys_info():

Qt WebEngine seems to be initialized from a plugin. Please set Qt::AA_ShareOpenGLContexts using QCoreApplication::setAttribute before constructing QGuiApplication.
Platform:      Linux-5.11.0-44-generic-x86_64-with-debian-bullseye-sid
Python:        3.7.9 (default, Aug 31 2020, 12:42:55)  [GCC 7.3.0]
Executable:    /home/katharina/snap/anaconda3/envs/py37cmp-gui/bin/python
CPU:           x86_64: 16 cores
Memory:        31.3 GB

mne:           0.20.7
numpy:         1.18.5 {blas=mkl_rt, lapack=mkl_rt}
scipy:         1.5.2
matplotlib:    3.2.2 {backend=Qt5Agg}

sklearn:       1.0.1
numba:         Not found
nibabel:       3.2.1
cupy:          Not found
pandas:        1.3.4
dipy:          1.1.1
mayavi:        Not found
pyvista:       Not found
vtk:           Not found

[welcome[bot]]

Hello! 👋 Thanks for opening your first issue here! ❤️ We will try to get back to you soon. 🚴🏽‍♂️

[agramfort]

@olafhauk can you have a look?

🙏 @Katharinski for the bug report

[drammock]

cross-ref to related forum post: https://mne.discourse.group/t/computing-the-resolution-matrix-with-loose-orientations/4148

[pettetmw]

Apologies for the long-winded post, but I hope it will help move the conversation forward. I've done a little digging around in the resolution matrix source code, and have put together the following notes. I doubt that everything is accurate or clear, so I welcome any feedback that helps refine the content into something useful as a specification for a workaround or patch. I've flagged a few points of uncertainty with "(?)". Also, at the end is a question about how to obtain the desired n_dipoles-by-n_dipoles shape for the resolution matrix in the free orientation case.

It appears to me that mne.minimum_norm.make_inverse_resolution_matrix() intends to return n_dipoles-by-n_dipoles only when both the forwardand inverse_operatorarguments have "fixed" (i.e. constrained to surface normal) orientation. If they are "free" orientation, the inverse is built from a VectorSourceEstimate, whose data property is "array of shape (n_dipoles, 3, _nsensors). Each dipole contains three vectors that denote the dipole strength in X, Y and Z directions over sensors." (Note that the algorithm spoofs the time dimension of the VectorSrouceEstimate object to represent the sensor dimension.) In this case, later in the process:

if invmat.ndim == 3:
    shape = invmat.shape
    invmat = invmat.reshape(shape[0] * shape[1], shape[2])

which should be (n_dipoles 3)-by-n_sensors, where the rows are the three x,y,z dimensions nested in the n_dipoles: [ x0, y0, z0, x1, y1, z1, ... etc ]. The same appears to be true(?) for the nesting of the columns in leadfield (=forward['so']['data']), which is n_sensors-by-(n_dipoles 3). So, the dot product,

invmat.dot(leadfield)

should be (n_dipoles x 3)-by-(n_dipoles x 3). This can be viewed as an n_dipoles-by-n_dipoles matrix of 3-by-3 matrices. Although a bit mind-bending, I think this might actually be logical when considering "free", or unconstrained, orientation in the context of the resolution matrix. We estimate CT/PS between any two unit-valued dipoles as the magnitude of the projection (?) of one onto the other caused by applying the inverse to the leadfield (superposed over all the sensors). When the dipole orientations are constrained, we obtain one projection coefficient for each dipole pair. (The projection is non-commutative(?), so the n_dipoles-by-n_dipoles matrix of pairwise coefficients is not symmetric(?)) On the other hand, when orientation is unconstrained, each 3-by-3 matrix is some kind of tensor-like thing(?), where each entry represents the pairwise projections between x, y, and z unit vectors at the two dipole locations. In other words, each matrix computes a 3-d dipole estimate at one location from a 3-d dipole generator at another location. Given a constraining set of dipole orientations, these matrices can determine the scalar magnitude of the projections between all the dipole pairs, as in the case of the "fixed" orientation forward and inverse above. (again non-commutative and non-symmetric(?))

So, here's the question I mentioned earlier: assuming the foregoing is roughly correct for free orientations, can we summarize the CT/PS for each dipole location pair by some scalar descriptor of each 3-by-3 tensor thing? Maybe just the matrix trace?

[larsoner]

should be (n_dipoles x 3)-by-(n_dipoles x 3). This can be viewed as an n_dipoles-by-n_dipoles matrix of 3-by-3 matrices. Although a bit mind-bending, I think this might actually be logical when considering "free", or unconstrained, orientation in the context of the resolution matrix.

Agreed this should be what is effectively computed in the loose and free orientation cases

can we summarize the CT/PS for each dipole location pair by some scalar descriptor of each 3-by-3 tensor thing? Maybe just the matrix trace?

My vote would be to add a combine_xyz='spectral' argument that mirrors compute_depth_prior. Even for our most compact recommended source space, which has ~8000 sources, a resulting non-3x3-normed matrix would end up being around shape (24000, 24000), requiring over 4GB of memory, compared to the collapsed version which "only" requires ~500 MB. I would see how far we can get with combine_xyz always being either 'spectral' or 'fro', and assume that getting the whole 3x3 is YAGNI for now.

[larsoner]

... also FYI I realized that this is very closely related to https://github.com/mne-tools/mne-python/issues/8737 and https://github.com/mne-tools/mne-python/issues/8620, and there is some progress in https://github.com/mne-tools/mne-python/pull/8639. @pettetmw maybe you could try #8620 and see if it works for your use case?

[larsoner]

... sorry I meant try #8639

[pettetmw]

So in this case, combine_xyz='spectral' would imply something like?:

res_mat = np.linalg.norm(np.einsum('svj,svk->vjk', M, G),  # vector dot prods
                           ord=2, axis=(1, 2))  # ord=2 spectral (largest s.v.)
                 # or maybe 'vsj,svk->vjk' ?

where: res_mat is p x p resolution matrix M is p x n x 3 inverse operator matrix (whitened) G is n x p x 3 leadfield gain matrix (whitened) p is number of dipoles n is number of sensors

[pettetmw]

As for memory use, for now, I'm only interested in CTPS between a few labelled regions, so I can pick dipoles from M and G to reduce memory/computation load.

[pettetmw]

8639 looks promising. How about if we hold position on einsumthread until I can give it a shot?

[larsoner]

Yeah I think getting #8639 into shape is the way to go

[Katharinski]

Hello!

Sorry for the long silence - I'm still on this! Thank you for the discussion, it really helped me to understand that a resolution matrix of (3ndipoles) x (3ndipoles) is the correct way to go when handling free orientations.

I tried the solution in pull request #8639, but I think it doesn't entirely solve the problem. Yes, in _get_psf_ctf, there is a bit that reduces the vector to a scalar (starting at line 150 in resolution_matrix.py), but it only does that in one direction and I still end up with a resolution matrix of ndipoles x (3*ndipoles). Unless I'm missing something?

The solution I am using right now is to compute the full matrix and then take the spectral norm of each 3x3 matrix, as suggested by @larsoner and @pettetmw It's quite slow and not very memory efficient, but I don't understand that einsum business well enough to make a better suggestion :grimacing: and it works for now.

So, thanks for your help! I'll keep an eye on this issue, maybe the combine_xyz='spectral' argument will find its way into a new version :grin:

[pettetmw]

My compliments, and thanks for proof of concept. So, did you do something like, e.g. from compute_depth_prior() :

        # n_pos = G.shape[1] // 3
        # The following is equivalent to this, but 4-10x faster
        # d = np.zeros(n_pos)
        # for k in range(n_pos):
        #     Gk = G[:, 3 * k:3 * (k + 1)]
        #     x = np.dot(Gk.T, Gk)
        #     d[k] = linalg.svdvals(x)[0]

... except using np.dot(Mk.T,Gk), with "M" and "G" defined something like my Feb 21 post? No need to share code, if that's not convenient. I'm just looking for a thumb's up that I'm on the right track.

[Katharinski]

Yes! I don't have the code at hand right now, would share it otherwise, but that's exactly right. Not using the einsum, but the spectral norm, yes.

[olafhauk]

Sorry that I didn't join this discussion before. I'm now looking into this as part of a coding spring at the CBU. But it looks like you are on top of things. It seems that the documentation wasn't updated: The resolution matrix should be "n_inversen_forward", where the "n" reflects all components of the source model, i.e. 3n_dipoles for free/loose orientations. Please let me know if you have any new suggestions. I'll start by checking the documentation and then see if anything else can be done better.