Série N°3 : Administration système -- Correction

Ce TP porte sur la résolution du problème inverse en EEG par régularisation. Il présente deux méthodes principales : la régularisation de Tikhonov (norme L2) et la régularisation TV-L1, appliquées respectivement via les algorithmes MNE et SISSY. L’objectif est d’implémenter ces algorithmes, d’étudier l’influence des paramètres de régularisation et de comparer différentes heuristiques pour leur choix. Pour réaliser ce TP, il faut disposer des données EEG simulées, des matrices de gain A et de l’opérateur variation T, ainsi que d’un environnement de calcul permettant la manipulation matricielle et la visualisation 3D.

Objectifs

• Comprendre et implémenter la régularisation de Tikhonov (MNE) pour le problème inverse en EEG.
• Comprendre et implémenter la régularisation TV-L1 (SISSY) via l’algorithme ADMM.
• Analyser l’influence des paramètres de régularisation sur la qualité de la solution inverse.
• Tester et comparer différentes heuristiques pour le choix du paramètre de régularisation.
• Visualiser les solutions obtenues et interpréter les résultats.

Prérequis et préparation

• Connaissances de base en optimisation et en algèbre linéaire.
• Notions sur les problèmes inverses et la régularisation.
• Environnement de calcul avec fonctions matricielles et visualisation 3D (ex. MATLAB ou Python avec bibliothèques appropriées).
• Accès aux matrices de gain A, aux données x, et à l’opérateur variation T (calculable via la fonction variation operator(mesh, 'face')).

Algorithme MNE : régularisation de Tikhonov

Le problème inverse régularisé de Tikhonov s’écrit :

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

La solution analytique est donnée par :

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

où I est la matrice identité. Cette formule permet d’obtenir directement la solution sans itération.

Étapes :

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

Un résultat correct montre une solution lisse qui équilibre la fidélité aux données et la régularisation.

Étude de l’algorithme MNE

Pour approfondir la compréhension :

1. Faire varier λ entre 0.1 et 1000. Observer comment la solution s’éloigne ou se rapproche de la configuration source originale.
2. Faire varier le rapport signal sur bruit (RSB) entre 0.1 et 10, puis λ entre 1 et 1000. Comparer les résultats pour comprendre l’impact du bruit sur le choix de λ.
3. Tester trois heuristiques pour choisir λ, avec RSB fixé à 1 et λ variant entre 1e-10 et 1e10 :
• Critère de la courbe en L (L-curve) : tracer ||x − As||2 en fonction de ||s||2 sur une double échelle logarithmique. Choisir λ au pli de la courbe.
• Principe de la discrepancy : choisir λ pour que 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 λ qui minimise la fonction GCV(λ) définie par :

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

Comparer les performances des trois heuristiques et conclure sur leur efficacité dans le contexte EEG.

Algorithme SISSY : régularisation TV-L1 via ADMM

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 variation (gradient sur la surface du maillage cortex). 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 associé à la norme L1, 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. 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érations (fixée à 60), mettre à jour s par :

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

Étapes :

1. Écrire la fonction SISSY(x, A, T, lambda, alpha) implémentant ADMM avec les mises à jour ci-dessus.
2. Tester pour RSB = 10, λ = 1, α = 0.1.
3. Visualiser la solution et comparer avec la configuration source originale.

Étude de l’algorithme SISSY

1. Faire varier λ entre 0.01 et 1000 et observer l’impact sur la solution inverse.
2. Fixer une valeur appropriée de λ, puis faire varier α de 0 à 1. Comparer les solutions obtenues à la configuration source.
3. Conclure sur l’influence respective des paramètres λ et α sur la qualité de la solution inverse.
4. Pour choisir λ, considérer la contrainte réelle que l’on souhaite imposer :

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

où ||·||_0 est la norme L0 (nombre d’éléments non nuls). Pour calculer cette norme, considérer que tous les éléments s_d tels que :

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

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

Résultats attendus

• Pour MNE, la solution doit être stable et lisse, avec un compromis entre fidélité aux données et régularisation contrôlé par λ.
• La variation de λ montre que une régularisation trop faible amplifie le bruit, tandis qu’une régularisation trop forte lisse excessivement la solution.
• Les heuristiques L-curve, discrepancy principle et GCV doivent permettre de choisir un λ pertinent, avec des performances variables selon le contexte.
• Pour SISSY, la solution est plus parcimonieuse et structurée, grâce à la régularisation TV-L1.
• Les paramètres λ et α influencent la sparsité et la structure de la solution. Un bon choix équilibre ces effets.
• L’heuristique basée sur la norme L0 permet de sélectionner λ en minimisant la complexité réelle de la solution.

Pièges courants

• Confondre les indices dans la définition de l’opérateur T, ce qui fausse la régularisation TV-L1.
• Ne pas utiliser la décomposition de Cholesky dans SISSY, ce qui rend le calcul de s très coûteux.
• Choisir un nombre d’itérations ADMM insuffisant (moins de 60), empêchant la convergence.
• Oublier d’ajuster le paramètre ρ (ici fixé à 1) dans ADMM, ce qui peut ralentir la convergence.
• Ne pas normaliser ou vérifier la puissance du bruit pour appliquer correctement le principe de discrepancy.
• Ne pas visualiser les solutions, ce qui empêche de détecter des artefacts ou des erreurs dans l’implémentation.