Back to the main page.

Bug 3447 - Implement regularization in reduced subspace as discussed at BIOMAG

Reported 2018-09-04 10:47:00 +0200
Modified 2019-04-01 08:57:59 +0200
Product: FieldTrip
Component: inverse
Version: unspecified
Hardware: PC
Operating System: Windows
Importance: P5 major
Assigned to: Jan-Mathijs Schoffelen
Depends on:
See also:

Vladimir Litvak - 2018-09-04 10:47:08 +0200

As discussed with Olaf Hauk and others at BIOMAG, the proper way to regularize for MaxFiltered data would be in reduced subspace after removing components with zero eigenvalues. I'm not sure how you want this implemented. If you get this working I can test on the phantom data and Krish Singh can test on his whole brain connectivity example.

Vladimir Litvak - 2018-09-13 19:20:29 +0200

Here is a function from Olaf to do combined matrix inversion with regularization when the rank is known. I could hack it into beamformer_lcmv etc. just by replacing the line with pinv and test. But you also have the whole subspace case there so you might prefer a different solution: function matreg = Tikhonov_rank_def(mat, rank, lambda) % Apply Tikhonov regularisation to rank-deficient matrix % mat: square matrix % rank: number of singular values to be considered in inversion % lambda: regularisation parameters % OH, Sep 2018 % SVD of input matrix [U,S,V] = svd(mat); % get singular values s = diag(S); % take only relevant values s2 = s(1:rank); % regularise eigenvalues with Tikhonov and invert s2 = s2 ./ (s2.^2 + lambda); % reconstitute regularised inverse matrix matreg = V(:,1:rank)*diag(s2)*U(:,1:rank)'; Test that the function produces the desired result >> mat = randn(3,3) mat = 2.7694 0.7254 -0.2050 -1.3499 -0.0631 -0.1241 3.0349 0.7147 1.4897 >> matreg = Tikhonov_rank_def(mat, 3, 0.1) matreg = 0.1480 -0.4172 -0.0038 0.5259 1.4010 0.1864 -0.5473 0.1265 0.5609 >> mat'*inv(mat*mat' + 0.1*eye(3)) ans = 0.1480 -0.4172 -0.0038 0.5259 1.4010 0.1864 -0.5473 0.1265 0.5609

Robert Oostenveld - 2018-09-14 09:07:25 +0200

I discussed at BIOMAG with Britta and Sarang (now CC) that we should have an "automatic cliff detection" to determine/estimate the numerical rank of the covariance. The cliff refers to the sudden drop in the singular value spectrum. There were some cases that we identified as resulting in extra cliffs. - Combined planar/magnetometer - Combined eeg/meg - Maxfiltered data - ica cleaned data - Concatenated data from different fif files that have been max filtered separately - average referenced eeg (or bipolar, or another silly reference scheme)

Jan-Mathijs Schoffelen - 2019-02-05 11:32:05 +0100

There is now a dedicated ft_inv function, which has replaced the call to pinv in the respective beamformer functions, for covariance inversion. It now supports various flavours of inversion: -truncated svd -> truncation based on 'kappa', or 'tolerance' parameter (tolerance is known in the cfg as 'tol', due to potential clashes with a 'tolerance' option that has a different meaning. -diagonally loaded (with lambda) pinv (original only option, still default), method is called 'tikhonov' -simple moorepenrose inversion inv(X'*X)*X' -simple invers inv() -truncated svd with replacement of low singular values (this is most akin to Olaf's implementation), method called 'winsorize' -> truncation based on 'kappa', or 'tolerance', all remaining singular values < lambda replaced by lambda. To be tested, documented and broadcasted to the community...