This is an automated email from the git hooks/post-receive script. yoh pushed a commit to tag 0.4 in repository python-mne.
commit 903a20c14ab0aa99a56217fc02cdd3c2333c9be3 Author: Alexandre Gramfort <[email protected]> Date: Wed Jul 18 09:43:38 2012 +0200 ENH : add better debiasing as suggested by Daniel --- mne/mixed_norm/debiasing.py | 116 +++++++++++++++++++++++++++++++++ mne/mixed_norm/optim.py | 9 ++- mne/mixed_norm/tests/test_debiasing.py | 22 +++++++ 3 files changed, 142 insertions(+), 5 deletions(-) diff --git a/mne/mixed_norm/debiasing.py b/mne/mixed_norm/debiasing.py new file mode 100755 index 0000000..e962a4e --- /dev/null +++ b/mne/mixed_norm/debiasing.py @@ -0,0 +1,116 @@ +# Authors: Daniel Strohmeier <[email protected]> +# Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +from math import sqrt +import numpy as np +from scipy import linalg + +from ..utils import check_random_state + + +def power_iteration_kron(A, C, max_iter=1000, tol=1e-3, random_state=0): + """Find the largest singular value for the matrix kron(C.T, A) + + It uses power iterations. + + Parameters + ---------- + A : array + An array + C : array + An array + max_iter : int + Maximum number of iterations + random_state : int | RandomState | None + Random state for random number generation + + Returns + ------- + L : float + largest singular value + """ + AS_size = C.shape[0] + rng = check_random_state(random_state) + B = rng.randn(AS_size, AS_size) + B /= linalg.norm(B, 'fro') + ATA = np.dot(A.T, A) + CCT = np.dot(C, C.T) + L0 = np.inf + for _ in range(max_iter): + Y = np.dot(np.dot(ATA, B), CCT) + L = linalg.norm(Y, 'fro') + + if abs(L - L0) < tol: + break + + B = Y / L + L0 = L + return L + + +def compute_bias(M, G, X, max_iter=1000, tol=1e-4, n_orient=1): + """Compute scaling to correct amplitude bias + + It solves the following optimization problem using FISTA: + + min 1/2 * (|| M - GDX ||fro)^2 + s.t. D >= 0 and D is a diagonal matrix + + Parameters + ---------- + M : array + measurement data + G : array + leadfield matrix + X : array + reconstructed time courses with amplitude bias + max_iter : int + Maximum number of iterations + tol : float + The tolerance on convergence + n_orient : int + The number of orientations (1 for fixed and 3 otherwise) + + Returns + ------- + D : array + Debiasing weights + """ + n_sources = X.shape[0] + + lipschitz_constant = 1.1 * power_iteration_kron(G, X) + + # initializations + D = np.ones(n_sources) + Y = np.ones(n_sources) + t = 1.0 + + for i in xrange(max_iter): + D0 = D + + # gradient step + R = M - np.dot(G * Y, X) + D = Y + np.sum(np.dot(G.T, R) * X, axis=1) / lipschitz_constant + # Equivalent but faster than: + # D = Y + np.diag(np.dot(np.dot(G.T, R), X.T)) / lipschitz_constant + + # prox ie projection on constraint + if n_orient != 1: # take care of orientations + # The scaling has to be the same for all orientations + D = np.mean(D.reshape(-1, n_orient), axis=1) + D = np.tile(D, [n_orient, 1]).T.ravel() + D = np.maximum(D, 0.0) + + t0 = t + t = 0.5 * (1.0 + sqrt(1.0 + 4.0 * t ** 2)) + Y.fill(0.0) + dt = (t0 - 1.0) / t + Y = D + dt * (D - D0) + if linalg.norm(D - D0, np.inf) < tol: + print "Debiasing converged after %d iterations" % i + break + else: + print "Debiasing did not converge" + return D diff --git a/mne/mixed_norm/optim.py b/mne/mixed_norm/optim.py index efe5ac3..147c114 100644 --- a/mne/mixed_norm/optim.py +++ b/mne/mixed_norm/optim.py @@ -6,6 +6,8 @@ from math import sqrt import numpy as np from scipy import linalg +from .debiasing import compute_bias + def groups_norm2(A, n_orient): """compute squared L2 norms of groups inplace""" @@ -295,10 +297,7 @@ def mixed_norm_solver(M, G, alpha, maxit=200, tol=1e-8, verbose=True, n_orient=n_orient) if (active_set.sum() > 0) and debias: - # Run ordinary least square on active set - GX = np.dot(G[:, active_set], X) - k, _, _, _ = linalg.lstsq(GX.reshape(-1, 1), M.reshape(-1, 1)) - X *= k - GX *= k + bias = compute_bias(M, G[:, active_set], X, n_orient=n_orient) + X *= bias[:, np.newaxis] return X, active_set, E diff --git a/mne/mixed_norm/tests/test_debiasing.py b/mne/mixed_norm/tests/test_debiasing.py new file mode 100755 index 0000000..6ecda6d --- /dev/null +++ b/mne/mixed_norm/tests/test_debiasing.py @@ -0,0 +1,22 @@ +# Authors: Daniel Strohmeier <[email protected]> +# Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +import numpy as np +from numpy.testing import assert_almost_equal + +from ..debiasing import compute_bias + + +def test_compute_debiasing(): + """Test source amplitude debiasing""" + rng = np.random.RandomState(42) + G = rng.randn(10, 4) + X = rng.randn(4, 20) + debias_true = np.arange(1, 5, dtype=np.float) + M = np.dot(G, X * debias_true[:, np.newaxis]) + debias = compute_bias(M, G, X, max_iter=10000, n_orient=1, tol=1e-7) + assert_almost_equal(debias, debias_true, decimal=5) + debias = compute_bias(M, G, X, max_iter=10000, n_orient=2, tol=1e-5) + assert_almost_equal(debias, [1.8, 1.8, 3.72, 3.72], decimal=2) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
