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

||x − As||2

2 + λ||s||2

2

s = AT(AAT + λI)−1x.

Remarque : utiliser le lemme d’inversion (BTR−1B + P−1)−1BTR−1 = PBT(BPBT + R)−1.

(3)

(4)

1

2.2 Algorithme SISSY

A partir du problème d’optimisation

min

s

||x − As||2

2 + λ(||z||1 + α||y||1)

s. t. z = Ts, y = s

(5)

montrer que les équations de mise à jour de l’algorithme ADMM pour les variables s, z, y et les multiplicateurs

de Lagrange u et v sont donnés par :

ATx + ρTTz(k) + TTu(k) + ρy(k) + v(k)(cid:17)

s(k+1) = (cid:0)ATA + ρ(TTT + I)(cid:1)−1 (cid:16)

(cid:19)

(cid:18)

z(k+1) = proxλ/ρ

y(k+1) = proxλα/ρ

(cid:16)

u(k+1) = u(k) + ρ

v(k+1) = v(k) + ρ

(cid:16)

(cid:18)

(cid:19)

u(k)

Ts(k+1) −

1

ρ

1

s(k+1) −

ρ

z(k+1) − Ts(k+1)(cid:17)

y(k+1) − s(k+1)(cid:17)

v(k)

(6)

(7)

(8)

(9)

(10)

où proxβ(y) = minx

est connue :

1

2 ||x − y||2

2 + β||x||1 est l’opérateur proximal associé avec la norme L1 et dont la solution

xi =







yi − β si yi > β

yi + β si yi < −β

0

sinon

(11)

et ρ est un paramètre de pénalité (utiliser ρ = 1).

3 Manipulations

3.1 Implémentation de l’algorithme MNE

Ecrire une fonction MNE(x,A,lambda) qui implémente l’algorithme MNE. Tester l’algorithme pour λ = 1.

Visualiser la solution obtenue à l’aide de la fonction trisurf.

3.2 Etude de l’algorithme MNE

1. Faire varier le parametre de régularisation (entre 0.1 et 1000) et observer l’influence de ce parametre

sur le résultat de l’algorithme en comparant la solution inverse obtenue avec la configuration de source

originale.

2. Faire varier le RSB (entre 0.1 et 10) et faire varier le paramètre de régularisation (entre 1 et 1000).

Comparer les résultats. Que pouvez-vous conclure sur le choix du paramètre de régularisation adéquat

en fonction du RSB ?

3. Dans la littérature, on trouve plusieurs heuristiques qui sont utilisées pour fixer le paramètre de

régularisation. Par la suite, nous allons considérer trois de ces critères :

L-curve criterion Cette heuristique consiste à tracer l’erreur de reconstruction ||x−As||2 en fonction

de la norme L2 de la solution ||s||2 pour différents paramètres de régularisation sur une double échelle

logarithmique. La courbe ainsi obtenue est souvent en forme de L. Le paramètre de régularisation

est alors choisi de maniere a correspondre au pli du L.

Discrepancy principle Cette heuristique est basée sur l’idée que l’erreur de reconstruction est due

au bruit. Le paramètre de régularisation est alors choisie tel que la puissance de l’erreur de recons-

truction correspond a peu pres a la puissance du bruit, c’est-a-dire ||x − As||2

2. La puissance

du bruit doit, en pratique, être estimée préalablement à partir des données (mais elle est connue

pour notre exemple de données simulées).

2 ≈ ||n||2

2

Generalized cross-validation Cette heuristique est basée sur l’idée qu’un bon choix du paramètre

de régularisation devrait permettre de bien prédire des observations exclues dans le calcul de la

solution inverse. Suivant le principe de validation croisée, le paramètre de régularisation est alors

choisi tel qu’il minimise la fonction suivante qui caractérise l’erreur de prédiction généralisée :

GCV(λ) =

||x − As||2

2

(trace{I − A(ATA + λI)−1AT})2 =

||x − As||2

2

(trace{I − AAT(AAT + λI)−1})2 .

Tester les trois heuristiques pour le choix du parametre de régularisation en fixant le RSB a 1 et

en faisant varier le paramètre de régularisation entre 1e-10 et 1e10. Conclure sur l’efficacité des trois

heuristiques dans le contexte du problème inverse en EEG.

3.3 Implémentation de l’algorithme SISSY

Ecrire une fonction SISSY(x,A,T,lambda,alpha) qui implémente l’algorithme ADMM pour résoudre le

probleme d’optimisation (5). La matrice T peut être calculé a l’aide de la fonction variation operator

(mesh,’face’). Pour diminuer la complexité numérique dans le calcul de la mise à jour du vecteur s,

une décomposition matricielle (décomposition de Cholesky) et un lemme d’inversion sont utilisés. De ce fait,

insérer les lignes suivantes au début du code (exécutées une fois au début de l’algorithme) :

P=sparse(rho(T.’T+speye(size(T,2))));

APi=A/P;

L=chol(eye(size(A,1))+APi*A.’,’lower’);

s1=A.’*x;

et inclure les lignes suivantes dans la boucle pour la mise a jour du vecteur s a chaque itération de l’ADMM :

b=s1+rho(T.’(z+u/rho)+y+v/rho);

s=P\b-APi’(L’\(L\(APib)));

Fixer le nombre d’itérations de l’ADMM à 60.

Tester l’algorithme pour RSB=10, λ = 1, α = 0.1.

3.4 Etude de l’algorithme SISSY

1. Faire varier le parametre de régularisation λ (entre 0.01 et 1000) et observer l’influence de ce parametre

sur le résultat de l’algorithme.

2. Choisir une valeur approprié pour le parametre λ et faire varier le parametre α de 0 à 1. Comparer la

solution inverse avec la configuration de sources originale.

3. Conclure sur l’influence des deux paramètres de régularisation sur la solution inverse et du choix

approprié de ces paramètres.

4. Pour choisir un parametre de régularisation approprié, trouver le parametre λ pour lequel la contrainte

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

que nous souhaitons réellement imposer (basée sur la norme L0 et non la norme L1) est minimale. Pour

calculer la norme L0, considérer que tous les éléments pour lesquels |sd| < 0.01 maxd∈{1,...,D}(|sd|) sont

zéros. Tester cette heuristique pour le choix de λ et comparer avec le résultat obtenu en utilisant le

discrepancy principle.

3