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

Ce TP aborde la résolution du problème inverse en EEG par régularisation, en comparant deux méthodes : 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, afin de comprendre l’impact des paramètres de régularisation sur la qualité de la reconstruction des sources corticales. Pour réaliser ce TP, il est nécessaire de disposer d’un environnement de calcul capable de manipuler des matrices et d’exécuter des algorithmes itératifs, ainsi que des données simulées ou réelles d’EEG et des maillages du cortex.

Objectifs

• Implémenter l’algorithme MNE basé sur la régularisation de Tikhonov (norme L2).
• Implémenter l’algorithme SISSY basé sur la régularisation TV-L1 et l’algorithme ADMM.
• Analyser l’influence des paramètres de régularisation sur 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 de sources originale.

Prérequis et préparation

• Connaissances en optimisation, notamment sur les problèmes régularisés et les normes L1 et L2.
• Maîtrise des algorithmes itératifs, en particulier ADMM (Alternating Direction Method of Multipliers).
• Environnement de calcul capable de manipuler des matrices creuses et d’exécuter des décompositions de Cholesky.
• Données EEG simulées ou réelles, ainsi que la matrice de leadfield A et l’opérateur de variation T calculé à partir du maillage cortical.
• Fonction de visualisation trisurf pour afficher les solutions sur le cortex.

Régularisation de Tikhonov : algorithme MNE

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

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, A la matrice de leadfield, x les données EEG, et I la matrice identité.

Étapes :

• Écrire une fonction MNE(x, A, lambda) qui calcule s selon la formule ci-dessus.
• Tester la fonction pour λ = 1.
• Visualiser la solution obtenue avec la fonction trisurf pour vérifier que la source estimée est plausible.

Un résultat correct montre une reconstruction cohérente des sources corticales avec un équilibre entre fidélité aux données et régularisation.

Étude de l’algorithme MNE

Pour comprendre l’impact du paramètre λ :

• Faire varier λ entre 0.1 et 1000 et observer l’effet sur la solution inverse comparée à la source originale.
• Faire varier le rapport signal sur bruit (RSB) entre 0.1 et 10, puis λ entre 1 et 1000, et comparer les résultats.
• Tester trois heuristiques pour choisir λ avec un RSB fixé à 1 et λ variant entre 1e-10 et 1e10 :
• Critère de la courbe en L : tracer ||x − As||2 en fonction de ||s||2 sur une échelle logarithmique double. Le pli de la courbe indique le λ optimal.
• Principe de la discrepancy : choisir λ pour que la puissance de l’erreur ||x − As||2^2 soit proche de la puissance du bruit ||n||2^2.
• Validation croisée généralisée (GCV) : minimiser la fonction GCV(λ) = ||x − As||2^2 / (trace{I − A (AᵀA + λI)⁻¹ Aᵀ})^2.

Ces tests permettent de conclure sur l’efficacité relative des heuristiques dans le contexte EEG.

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

avec T l’opérateur gradient sur la surface du maillage cortical, λ et α les paramètres de régularisation, et les variables auxiliaires z, y.

L’algorithme ADMM met à jour les variables s, z, y et les multiplicateurs de Lagrange u, v selon :

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))

où prox_β(y) est l’opérateur proximal de la norme L1 défini par :

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

Le paramètre ρ est fixé à 1.

Implémentation de l’algorithme SISSY

Pour réduire la complexité numérique lors de la mise à jour de s, utiliser la décomposition de Cholesky et un lemme d’inversion :

P = sparse(ρ * (Tᵀ T + I))
APi = A / P
L = chol(I + APi * Aᵀ, 'lower')
s1 = Aᵀ x

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

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

Étapes :

• Écrire la fonction SISSY(x, A, T, lambda, alpha) implémentant ADMM selon les formules ci-dessus.
• Tester l’algorithme avec RSB = 10, λ = 1, α = 0.1.

Une solution correcte doit montrer une reconstruction plus parcimonieuse et structurée des sources corticales que MNE.

Étude de l’algorithme SISSY

Pour analyser l’impact des paramètres :

• Faire varier λ entre 0.01 et 1000 et observer l’effet sur la solution inverse.
• Fixer une valeur appropriée de λ, puis faire varier α de 0 à 1 et comparer les solutions obtenues à la configuration originale.
• Conclure sur l’influence conjointe de λ et α sur la qualité et la parcimonie de la solution.
• Pour choisir λ, chercher la valeur minimisant la contrainte réelle basée sur la norme L0 :

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

où un élément s_d est considéré nul si |s_d| < 0.01 max_d |s_d|.

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 montrer un compromis entre fidélité aux données et régularisation selon λ.
• La courbe en L doit apparaître lors du tracé de ||x − As||2 versus ||s||2, permettant d’identifier un λ pertinent.
• Les heuristiques de choix de λ doivent fournir des valeurs proches, avec des différences selon le bruit et la configuration.
• Pour SISSY, la solution doit être plus parcimonieuse, avec des sources localisées et moins bruitées.
• L’influence de α doit moduler la sparsité et la structure spatiale des sources.
• La minimisation de la norme L0 doit conduire à un choix de λ cohérent avec une bonne reconstruction.

Pièges courants

• Confondre les paramètres λ et α et leur rôle respectif dans la régularisation.
• Ne pas normaliser ou mal estimer le bruit, ce qui fausse le choix du paramètre de régularisation.
• Oublier d’utiliser la décomposition de Cholesky pour SISSY, ce qui rend le calcul de s très coûteux.
• Ne pas vérifier la convergence de l’algorithme ADMM (fixer un nombre d’itérations suffisant).
• Interpréter incorrectement l’opérateur proximal, notamment les seuils appliqués.
• Visualiser les résultats sans comparer à la configuration originale, ce qui empêche d’évaluer la qualité réelle de la reconstruction.