Administration système UNIX - Correction

Ce TP porte sur la résolution du problème inverse en EEG par régularisation, en mettant en œuvre deux méthodes : la régularisation de Tikhonov (norme L2) et la régularisation TV-L1. Il permet d’apprendre à coder et tester 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 d’un environnement de calcul capable d’exécuter du code matriciel (par exemple MATLAB ou Python avec les bibliothèques adéquates) et connaître les bases de l’optimisation et du traitement du signal.

Objectifs

• Implémenter l’algorithme MNE basé sur la régularisation de Tikhonov (norme L2).
• Étudier l’influence du paramètre de régularisation sur la solution inverse.
• Tester différentes heuristiques pour le choix du paramètre de régularisation.
• Implémenter l’algorithme SISSY utilisant la régularisation TV-L1 avec ADMM.
• Analyser l’impact des paramètres λ et α dans l’algorithme SISSY.
• Comparer les résultats obtenus avec les configurations de sources originales.

Prérequis et préparation

• Connaissances en optimisation convexe et en algèbre linéaire.
• Notions sur les problèmes inverses et la régularisation.
• Environnement de calcul capable de manipuler des matrices creuses et d’exécuter des algorithmes itératifs (MATLAB, Python).
• Fonctionnalités pour visualiser des surfaces 3D (fonction trisurf ou équivalent).
• Compréhension de la méthode ADMM et des opérateurs proximaux.

Introduction au problème inverse et régularisation

Le problème inverse en EEG consiste à estimer la source s à partir des données mesurées x, en résolvant l’optimisation :

min_s ||x − As||2^2 + λ f(s)

Le premier terme garantit la fidélité aux données, le second, appelé terme de régularisation, impose des contraintes sur la solution. Le paramètre λ équilibre ces deux termes.

Deux types de régularisation sont étudiés :

• Régularisation de Tikhonov (norme L2) : f(s) = ||s||2^2.
• Régularisation TV-L1 : f(s) = ||Ts||1 + α||s||1, où T est un opérateur gradient sur la surface du cortex.

L’opérateur T est défini par :

Te,d =
  1   si d = de,1
 -1   si d = de,2
  0   sinon

où de,1 et de,2 sont les indices des dipôles partageant l’arête e.

Algorithme MNE : solution analytique

Pour la régularisation L2, la solution analytique du problème :

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

est donnée par :

s = Aᵀ (A Aᵀ + λ I)^(-1) x

Cette formule est obtenue en utilisant un lemme d’inversion matricielle.

Algorithme SISSY : optimisation avec ADMM

Le problème d’optimisation avec régularisation TV-L1 s’écrit :

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

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

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

z^(k+1) = prox_(λ/ρ)(Ts^(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)_i =
  y_i − β   si y_i > β
  y_i + β   si y_i < −β
  0         sinon

Le paramètre ρ est fixé à 1.

Implémentation de l’algorithme MNE

Écrire une fonction MNE(x, A, lambda) qui calcule la solution analytique :

function s = MNE(x, A, lambda)
  s = A' * ((A * A' + lambda * eye(size(A,1))) \ x);
end

Tester cette fonction avec λ = 1. Visualiser la solution avec la fonction trisurf pour vérifier la reconstruction.

Étude de l’algorithme MNE

1. Faire varier le paramètre λ entre 0.1 et 1000. Observer comment la solution s’adapte, en comparant avec la configuration originale des sources.
2. Faire varier le rapport signal sur bruit (RSB) entre 0.1 et 10, puis λ entre 1 et 1000. Comparer les résultats et conclure sur le choix optimal de λ selon le RSB.
3. Tester trois heuristiques pour choisir λ avec RSB fixé à 1 et λ variant de 1e-10 à 1e10 :
• Critère de la courbe en L : tracer ||x − As||2 en fonction de ||s||2 sur une échelle logarithmique double. Choisir λ au pli de la courbe.
• Principe de la discrepancy : choisir λ tel que ||x − As||2^2 ≈ 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)^(-1) Aᵀ})^2.
4. Conclure sur l’efficacité relative de ces heuristiques dans le contexte EEG.

Implémentation de l’algorithme SISSY

Écrire une fonction SISSY(x, A, T, lambda, alpha) implémentant ADMM selon les mises à jour précédentes. Pour optimiser le calcul de s, insérer au début :

P = sparse(rho * (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 (60 itérations), mettre à jour s par :

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

Tester avec RSB = 10, λ = 1, α = 0.1.

É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 la solution avec la configuration originale.
3. Conclure sur l’influence des paramètres λ et α sur la solution et sur le choix optimal.
4. Pour choisir λ, minimiser la contrainte f(s) = ||Ts||0 + α ||s||0, où la norme L0 est approximée en considérant comme nuls tous les éléments s_d tels que |s_d| < 0.01 max_d |s_d|.
5. Tester cette heuristique et comparer avec le résultat obtenu par le principe de discrepancy.

Résultats attendus

• L’algorithme MNE doit fournir une solution stable et analytique pour tout λ, avec un compromis visible entre fidélité aux données et régularisation.
• La variation de λ modifie la lissage ou la sparsité de la solution : λ faible favorise la fidélité, λ élevé favorise la régularité ou la sparsité.
• Les heuristiques pour le choix de λ doivent converger vers des valeurs proches, avec des différences selon le niveau de bruit.
• L’algorithme SISSY doit converger en 60 itérations, produisant des solutions plus structurées et sparses grâce à la régularisation TV-L1.
• La variation de α influence la balance entre la pénalisation du gradient (||Ts||1) et la pénalisation directe (||s||1), modifiant la structure de la solution.
• L’heuristique basée sur la norme L0 approximée doit permettre de choisir un λ qui minimise la complexité effective de la solution.

Pièges courants

• Confondre les indices dans l’opérateur T, ce qui fausse la définition du gradient sur le maillage.
• Ne pas utiliser la décomposition de Cholesky et le lemme d’inversion dans SISSY, ce qui rend le calcul de s trop coûteux.
• Choisir un nombre d’itérations ADMM insuffisant, empêchant la convergence de SISSY.
• Ignorer la normalisation ou l’échelle des données, rendant les comparaisons entre solutions difficiles.
• Mal interpréter les heuristiques de choix de λ, notamment le critère de la courbe en L qui nécessite une échelle logarithmique et une lecture attentive du pli.
• Ne pas estimer correctement la puissance du bruit pour le principe de discrepancy, ce qui fausse le choix de λ.
• Omettre la mise à jour correcte des multiplicateurs de Lagrange u et v dans ADMM, ce qui empêche la convergence.