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

Ce TP explore la résolution du problème inverse en électroencéphalographie (EEG) par régularisation. Il met en œuvre deux méthodes principales : la régularisation de Tikhonov (norme L2) via l'algorithme MNE, et la régularisation TV-L1 via l'algorithme SISSY. Le TP permet de comprendre l'impact des paramètres de régularisation et d'expérimenter différentes heuristiques pour leur choix. Pour réaliser ce TP, il est nécessaire de disposer des données EEG simulées, des matrices de gain A et de l'opérateur gradient T, ainsi que d'un environnement de calcul capable d'exécuter des algorithmes matriciels et itératifs (par exemple MATLAB ou Python avec bibliothèques scientifiques).

Objectifs

• Comprendre et implémenter la régularisation de Tikhonov pour le problème inverse EEG.
• Mettre en œuvre et analyser l'algorithme ADMM pour la régularisation TV-L1 (SISSY).
• Étudier l'influence des paramètres de régularisation sur la qualité de la solution inverse.
• Tester différentes heuristiques pour le choix du paramètre de régularisation.
• Visualiser et comparer les solutions obtenues avec la configuration source originale.

Prérequis et préparation

• Connaissances de base en optimisation convexe et algorithmes itératifs.
• Maîtrise des notions de normes L1, L2, et opérateurs linéaires sur maillages.
• Environnement de calcul avec fonctions matricielles, décomposition de Cholesky et visualisation 3D (ex. MATLAB ou Python).
• Accès aux matrices A (matrice de gain EEG), T (opérateur gradient sur le maillage cortical), et aux données x (mesures EEG).

Régularisation de Tikhonov : algorithme MNE

Le problème inverse EEG est formulé comme :

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

La solution analytique s'écrit :

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

où λ est le paramètre de régularisation. Cette étape permet d'obtenir une solution stable en équilibrant la fidélité aux données et la norme de la solution.

À faire :

• Écrire une fonction MNE(x, A, lambda) qui calcule s selon la formule ci-dessus.
• Tester cette fonction avec λ = 1.
• Visualiser la solution s avec la fonction trisurf pour observer la distribution des sources sur le cortex.

Un résultat correct montre une solution lisse et stable, sans sur-ajustement aux données bruitées.

Étude de l'algorithme MNE

Cette étape vise à comprendre l'impact du paramètre λ et du rapport signal sur bruit (RSB) sur la solution inverse.

• Faire varier λ entre 0.1 et 1000 et observer les changements dans la solution inverse. Une faible λ favorise la fidélité aux données mais peut amplifier le bruit, tandis qu'une forte λ produit une solution plus lisse mais moins précise.
• Faire varier le RSB entre 0.1 et 10, puis pour chaque RSB, faire varier λ entre 1 et 1000. Comparer les résultats pour comprendre comment le bruit influence le choix optimal de λ.
• Tester trois heuristiques pour choisir λ :
• Critère de la courbe en L : tracer ||x − As||2 en fonction de ||s||2 sur une échelle logarithmique double. Le paramètre λ est choisi au "pli" de la courbe.
• Principe de la discrepancy : choisir λ tel que la puissance de l'erreur de reconstruction ||x − As||2^2 soit proche de la puissance du bruit ||n||2^2.
• Validation croisée généralisée (GCV) : choisir λ minimisant la fonction :

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

• Tester ces heuristiques avec RSB fixé à 1 et λ variant entre 1e-10 et 1e10.

Un bon choix de λ équilibre la qualité de la reconstruction et la complexité de la solution.

Régularisation TV-L1 : algorithme SISSY

Le problème d'optimisation est :

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

où T est l'opérateur gradient sur le maillage cortical, λ et α sont des paramètres de régularisation, et z, y sont des variables auxiliaires. L'algorithme ADMM est utilisé pour résoudre ce problème itérativement.

Les mises à jour des variables s, z, y, et des multiplicateurs de Lagrange u, v sont :

s^(k+1) = (AᵀA + ρ (TᵀT + I))⁻¹ (Aᵀ x + ρ Tᵀ (z^(k) − 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) + ρ (T s^(k+1) − z^(k+1))

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

L'opérateur proximal associé à la norme L1 est défini par :

prox_β(y)_i = 
  y_i − β si y_i > β
  y_i + β si y_i < −β
  0 sinon

avec ρ = 1.

Implémentation de l'algorithme SISSY

Pour réduire la complexité numérique lors de la mise à jour de s, on utilise une décomposition de Cholesky et un lemme d'inversion :

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érations ADMM (fixée à 60 itérations), la mise à jour de s est :

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.
• Tester avec RSB = 10, λ = 1, α = 0.1.

Une solution correcte montre une reconstruction plus parcimonieuse et structurée que MNE, grâce à la régularisation TV-L1.

Étude de l'algorithme SISSY

• Faire varier λ entre 0.01 et 1000 et observer son influence sur la solution inverse.
• Fixer une valeur appropriée de λ, puis faire varier α de 0 à 1. Comparer les solutions obtenues avec la configuration source originale.
• Conclure sur l'impact des paramètres λ et α sur la qualité et la structure de la solution inverse.
• Pour choisir λ, minimiser la contrainte :

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

où la norme L0 est approximée en considérant que tous les éléments |s_d| < 0.01 max_d(|s_d|) sont nuls.

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

Résultats attendus

• Pour MNE, la solution analytique doit être stable et lisse, avec une dépendance claire au paramètre λ.
• Les heuristiques de choix de λ doivent permettre d'identifier un compromis pertinent entre fidélité et régularisation.
• Pour SISSY, la solution doit être plus parcimonieuse et structurée, reflétant la régularisation TV-L1.
• L'algorithme ADMM doit converger en environ 60 itérations.
• La minimisation de la norme L0 approximée doit fournir un paramètre λ cohérent avec la qualité de la solution.

Pièges courants

• Confondre les indices dans l'opérateur T, ce qui fausse le calcul du gradient sur le maillage.
• Ne pas utiliser la décomposition de Cholesky pour SISSY, ce qui rend le calcul de s trop coûteux ou instable.
• Choisir un paramètre λ trop faible dans MNE, entraînant une solution bruitée et instable.
• Ne pas normaliser ou estimer correctement la puissance du bruit pour le principe de discrepancy.
• Oublier d'initialiser correctement les variables z, y, u, v dans l'algorithme ADMM.
• Ne pas vérifier la convergence de l'algorithme ADMM, ce qui peut conduire à des résultats erronés.