Problème inverse d'EEG : TP2 - Régularisation de Tikhonov et régularisation TV-L1

Ce TP porte sur la résolution du problème inverse en EEG à l’aide de deux méthodes de régularisation : la régularisation de Tikhonov (norme L2) et la régularisation TV-L1. Il permet d’implémenter et d’étudier ces deux algorithmes, MNE et SISSY, pour localiser des sources distribuées à partir de données EEG. Le travail nécessite des connaissances en optimisation, en algèbre linéaire, ainsi qu’une maîtrise de la programmation numérique (par exemple en MATLAB ou Python).

Objectifs

• Comprendre et implémenter la régularisation de Tikhonov (norme L2) pour le problème inverse EEG.
• Comprendre et implémenter la régularisation TV-L1 via l’algorithme ADMM (algorithme SISSY).
• Étudier l’influence des paramètres de régularisation sur la qualité de la solution inverse.
• Appliquer différentes heuristiques pour le choix du paramètre de régularisation.
• Visualiser et comparer les solutions obtenues avec la configuration de sources originale.

Prérequis et installation

• Connaissances en optimisation convexe et méthodes itératives (notamment ADMM).
• Maîtrise des notions de normes L1, L2, et TV (variation totale).
• Compétences en programmation numérique (MATLAB ou équivalent).
• Accès à une fonction de visualisation 3D telle que trisurf pour afficher les résultats.
• Connaissance du modèle EEG et des matrices associées (matrice de gain A, opérateur gradient T).

Régularisation de Tikhonov : algorithme MNE

L’objectif est de résoudre le problème inverse suivant :

min_s ||x − A s||2^2 + λ ||s||2^2

où x est le signal mesuré, A la matrice de gain, s la source inconnue, et λ le paramètre de régularisation.

La solution analytique est donnée par :

s = Aᵀ (A Aᵀ + λ I)⁻¹ x

avec I la matrice identité. Cette formule est obtenue en appliquant un lemme d’inversion spécifique.

À faire :

• Écrire une fonction MNE(x, A, lambda) qui calcule s selon la formule ci-dessus.
• Tester la fonction avec λ = 1.
• Visualiser la solution s avec la fonction trisurf pour vérifier la localisation des sources.

Un résultat correct doit montrer une reconstruction cohérente des sources sur le cortex modélisé.

Étude de l’algorithme MNE

Pour approfondir la compréhension :

1. Faire varier λ entre 0.1 et 1000. Observer comment la solution s change, notamment en termes de lissage et de fidélité aux données.
2. Faire varier le rapport signal sur bruit (RSB) entre 0.1 et 10, puis refaire varier λ entre 1 et 1000. Comparer les résultats pour comprendre l’impact du bruit sur le choix de λ.
3. Tester trois heuristiques pour choisir λ :
• Critère de la courbe en L : tracer sur une double échelle logarithmique l’erreur de reconstruction ||x − A s||2 en fonction de la norme ||s||2 pour différents λ. Le pli de la courbe indique le choix optimal.
• Principe de la discrepancy : choisir λ tel que la puissance de l’erreur de reconstruction soit proche de la puissance du bruit, c’est-à-dire ||x − A s||2^2 ≈ ||n||2^2.
• Validation croisée généralisée (GCV) : choisir λ minimisant la fonction :

GCV(λ) = ||x − A s||2^2 / (trace{I − A (Aᵀ A + λ I)⁻¹ Aᵀ})^2

4. Tester ces heuristiques avec RSB = 1 et λ variant de 1e-10 à 1e10. Comparer leur efficacité dans le contexte EEG.

Régularisation TV-L1 : algorithme SISSY

Le problème d’optimisation devient :

min_s ||x − A s||2^2 + λ (||z||1 + α ||y||1)
s.t. z = T s, y = s

où T est un opérateur gradient sur la surface du maillage cortical, z et y sont des variables auxiliaires, et α un paramètre supplémentaire.

L’algorithme ADMM est utilisé pour résoudre ce problème, avec les mises à jour :

s^(k+1) = (Aᵀ A + ρ (Tᵀ T + I))⁻¹ [Aᵀ x + ρ Tᵀ z^(k) + Tᵀ u^(k) + ρ y^(k) + v^(k)]

z^(k+1) = prox_(λ/ρ)(T s^(k+1) − u^(k)/ρ)

y^(k+1) = prox_(λα/ρ)(s^(k+1) − v^(k)/ρ)

u^(k+1) = u^(k) + ρ (z^(k+1) − T s^(k+1))

v^(k+1) = v^(k) + ρ (y^(k+1) − s^(k+1))

L’opérateur proximal prox_β(y) associé à la norme L1 est défini par :

prox_β(y) = argmin_x (1/2) ||x − y||2^2 + β ||x||1

avec la solution élémentaire :

x_i = 
  y_i − β   si y_i > β
  y_i + β   si y_i < −β
  0         sinon

Le paramètre de pénalité ρ est fixé à 1.

Implémentation de l’algorithme SISSY

Pour réduire la complexité numérique lors de la mise à jour de s, une décomposition de Cholesky et un lemme d’inversion sont utilisés. Au début de l’algorithme, insérer :

P = sparse(ρ * (T.' * T + speye(size(T, 2))));
APi = A / P;
L = chol(eye(size(A,1)) + APi * A.', 'lower');
s1 = A.' * x;

Dans la boucle d’itération (fixée à 60 itérations), mettre à jour s avec :

b = s1 + ρ * (T.' * (z + u / ρ) + y + v / ρ);
s = P \ b - APi.' * (L.' \ (L \ (APi * b)) );

À faire :

• Écrire une fonction SISSY(x, A, T, lambda, alpha) implémentant cet algorithme ADMM.
• Tester avec RSB = 10, λ = 1 et α = 0.1.

Une solution correcte doit montrer une localisation précise des sources avec un effet de régularisation favorisant la structure et la parcimonie.

Étude de l’algorithme SISSY

1. Faire varier λ entre 0.01 et 1000. Observer l’impact sur la solution inverse.
2. Fixer une valeur appropriée de λ, puis faire varier α de 0 à 1. Comparer les solutions obtenues avec la configuration originale des sources.
3. Conclure sur l’influence respective de λ et α dans la qualité de la solution.
4. Pour choisir un paramètre λ adapté, minimiser la contrainte :

f(s) = ||T s||_0 + α ||s||_0

où la norme L0 est approximée en considérant que tous les éléments s_d tels que :

|s_d| < 0.01 max_{d} |s_d|

sont nuls. Tester cette heuristique et comparer le résultat avec celui obtenu par le principe de discrepancy.

Résultats attendus

• Pour l’algorithme MNE, la solution analytique doit être stable et montrer un compromis entre fidélité aux données et lissage selon λ.
• Les heuristiques de choix de λ doivent permettre d’identifier un paramètre optimal, visible notamment par la courbe en L, la correspondance avec la puissance du bruit, ou la minimisation de la fonction GCV.
• Pour l’algorithme SISSY, la solution doit être plus parcimonieuse et structurée, avec une influence visible des paramètres λ et α sur la sparsité et la régularité.
• La minimisation de la contrainte basée sur la norme L0 doit permettre de choisir un λ qui favorise une solution plus proche de la réalité physique des sources.

Pièges courants

• Ne pas utiliser la bonne formule analytique pour MNE, notamment en inversant mal les matrices.
• Confondre les rôles des paramètres λ et α dans SISSY, ce qui peut conduire à une régularisation inadéquate.
• Ne pas normaliser ou estimer correctement la puissance du bruit, rendant le principe de discrepancy inefficace.
• Oublier de fixer un nombre suffisant d’itérations dans ADMM (au moins 60) pour assurer la convergence.
• Mal implémenter l’opérateur proximal, ce qui fausse la mise à jour des variables z et y.
• Ne pas visualiser les résultats avec une fonction adaptée, rendant l’interprétation difficile.