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

d’optimisation (1). Pour la plupart des termes de r´egularisation, comme, par exemple, pour la norme TV-L1,

ceci n´ecessite de recourir `a un algorithme it´eratif. Par contre, pour la norme L2, une solution analytique peut

ˆetre d´eduite.

2.1 Algorithme MNE

Montrer que la solution analytique au probl`eme de minimisation

peut s’´ecrire sous la forme :

min

s

||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`eme d’optimisation

min

s

||x − As||2

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

s. t. z = Ts, y = s

(5)

montrer que les ´equations de mise `a jour de l’algorithme ADMM pour les variables s, z, y et les multiplicateurs

de Lagrange u et v sont donn´es 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`u proxβ(y) = minx

est connue :

1

2 ||x − y||2

2 + β||x||1 est l’op´erateur proximal associ´e avec la norme L1 et dont la solution

xi =







yi − β si yi > β

yi + β si yi < −β

0

sinon

(11)

et ρ est un param`etre de p´enalit´e (utiliser ρ = 1).

3 Manipulations

3.1 Impl´ementation de l’algorithme MNE

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

Visualiser la solution obtenue `a l’aide de la fonction trisurf.

3.2 Etude de l’algorithme MNE

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

sur le r´esultat 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`etre de r´egularisation (entre 1 et 1000).

Comparer les r´esultats. Que pouvez-vous conclure sur le choix du param`etre de r´egularisation ad´equat

en fonction du RSB ?

3. Dans la litt´erature, on trouve plusieurs heuristiques qui sont utilis´ees pour fixer le param`etre de

r´egularisation. Par la suite, nous allons consid´erer trois de ces crit`eres :

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

de la norme L2 de la solution ||s||2 pour diff´erents param`etres de r´egularisation sur une double ´echelle

logarithmique. La courbe ainsi obtenue est souvent en forme de L. Le param`etre de r´egularisation

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

Discrepancy principle Cette heuristique est bas´ee sur l’id´ee que l’erreur de reconstruction est due

au bruit. Le param`etre de r´egularisation 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, ˆetre estim´ee pr´ealablement `a partir des donn´ees (mais elle est connue

pour notre exemple de donn´ees simul´ees).

2 ≈ ||n||2

2

Generalized cross-validation Cette heuristique est bas´ee sur l’id´ee qu’un bon choix du param`etre

de r´egularisation devrait permettre de bien pr´edire des observations exclues dans le calcul de la

solution inverse. Suivant le principe de validation crois´ee, le param`etre de r´egularisation est alors

choisi tel qu’il minimise la fonction suivante qui caract´erise l’erreur de pr´ediction g´en´eralis´ee :

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´egularisation en fixant le RSB a 1 et

en faisant varier le param`etre de r´egularisation entre 1e-10 et 1e10. Conclure sur l’efficacit´e des trois

heuristiques dans le contexte du probl`eme inverse en EEG.

3.3 Impl´ementation de l’algorithme SISSY

Ecrire une fonction SISSY(x,A,T,lambda,alpha) qui impl´emente l’algorithme ADMM pour r´esoudre le

probleme d’optimisation (5). La matrice T peut ˆetre calcul´e a l’aide de la fonction variation operator

(mesh,’face’). Pour diminuer la complexit´e num´erique dans le calcul de la mise `a jour du vecteur s,

une d´ecomposition matricielle (d´ecomposition de Cholesky) et un lemme d’inversion sont utilis´es. De ce fait,

ins´erer les lignes suivantes au d´ebut du code (ex´ecut´ees une fois au d´ebut 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´eration 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´erations de l’ADMM `a 60.

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

3.4 Etude de l’algorithme SISSY

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

sur le r´esultat de l’algorithme.

2. Choisir une valeur appropri´e pour le parametre λ et faire varier le parametre α de 0 `a 1. Comparer la

solution inverse avec la configuration de sources originale.

3. Conclure sur l’influence des deux param`etres de r´egularisation sur la solution inverse et du choix

appropri´e de ces param`etres.

4. Pour choisir un parametre de r´egularisation appropri´e, trouver le parametre λ pour lequel la contrainte

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

que nous souhaitons r´eellement imposer (bas´ee sur la norme L0 et non la norme L1) est minimale. Pour

calculer la norme L0, consid´erer que tous les ´el´ements pour lesquels |sd| < 0.01 maxd∈{1,...,D}(|sd|) sont

z´eros. Tester cette heuristique pour le choix de λ et comparer avec le r´esultat obtenu en utilisant le

discrepancy principle.

3