[PDF] [PDF] Mod`eles convexes et algorithmes doptimisation en imagerie

View - Download [PDF] Mod`eles convexes et algorithmes doptimisation en imagerie

Previous PDF

Next PDF

Modeles convexes et algorithmes d'optimisation

en imagerie.

Pierre Weiss.

October 25, 2011

III.2/ Theorie de la complexite en

optimisation convexe.

Applications a l'imagerie.

Algorithmes d'optimisation... Deuxieme partie.

Jusqu'a present, on n'est pas capable de resoudre les problemes annonces au depart !

Les descentes de sous-gradient sont trop lentes.

Les schemas acceleres ne s'adaptent pas aux

non-dierentiabilites.

Conclusion :

Pour obtenir des algorithmes plus ecaces, il faut : Se concentrer sur des classes de fonctions plus restreintes que les fonctions convexes quelconques.

Utiliser d'autres proprietes que les gradients et

sous-gradients (op erateurspro ximaux,dualit e). Quelques algorithmes d'optimisation non lisse ecaces. Deux mots sur les operateurs proximaux (Moreau 1960). Algorithmes primaux quand une partie du probleme est dierentiable. Algorithmes duaux quand une partie du probleme est fortement convexe.

Algorithmes primaux-duaux dans le cas general.

Operateurs proximaux.

[Combettes, Wajs 03], Signal recovery by proximal forward-backward splitting. SIAM MMS Denition :Soitf:Rn!R[f+1gun fonction convexe s.c.i.

On appelle

op erateurpro ximal ou r esolvante de fl'operateur : prox f(x0) = argmin x2Rnf(x) +12 kxx0k2:

Notes :

Cet operateur est bien deni carf(x) +12

kxx0k2est strictement convexe et coercive. Il est parfois note (Id +@f)1(x0) car les conditions d'optimalite menent a

02@f(proxf(x0)) + proxf(x0)x0:

Operateurs proximaux.

Intuition :

C'est une generalisation de la projection.

Sif(x) =0 six2X

+1sinon, alors : prox f(x0) = argmin x2X12 kxx0k2(1)

X(x0):(2)

Proprietes elementaires :

L'operateur proximal est non expansif :

8(x;y)2RnRn;kproxf(x)proxf(y)k kxyk:

Il caracterise les minimiseurs :

x2Argmin x2Rnf(x),02@f(x),proxf(x) =x(point xe):

Une classe de problemes non dierentiable :

L'operateur proximal est central pour la resolution - a l'aide de methodes de premier ordre - des problemes non dierentiables. Nous allons montrer son utilisation pour la resolution de quelques classes de probemes :

1.Somme d'une fonction convexe dierentiable et d'une

fonction convexe (Primal).

2.Somme d'une fonction convexe et d'une fonction fortement

convexe (Dual).

3.Autres problemes de type point selle.

Une classe de problemes non dierentiable Auslender70:

On considere desormais et jusqu'a

nouv elordre le probl eme suivant : minx2Rnf1(x) +f2(x) f1est convexe, dierentiable, a gradientL-Lispchitz. f2est convexe, s.c.i telle que proxf2peut ^etre calcule.

Propriete :Les minimiseursxsont caracterises par

x = proxf2(xrf1(x)):

Preuve :

02 rf1(x) +@f2(x)

,xrf1(x)2x+@f2(x) ,x= proxf2(xrf1(x)): Somme d'une fonction dierentiable et d'une fonction convexe.. Le gradient proximal :la caracterisation de point-xe precedente donne envie d'utiliser le schema : x k+1= proxf2(xkrf1(xk))

Exemple (compressive sensing) :

min x2Rn12 kAxbk2+kxk1 x k+1=shrink(xkAT(Axb)) [I. Daubechie etal 03] An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, CPAM.

Convergence lineaire :Si=min(ATA)>0, alors le

schema precedent converge lineairement en choisissant=2+L. Somme d'une fonction dierentiable et d'une fonction convexe..

Preuve :

kxk+1xk=kproxf2(xkAT(Axkb))xk =kproxf2(xkAT(Axkb)) proxf2(xAT(Axb))k k(IATA)(xkx)k(non expansif) jjj(IATA)jjj k(xkx)k avec :jjj(IATA)jjj= max(j1j;j1Lj).

Enn, on a

min max(1;1L) =2+L: Que se passe-t-il dans le cas non fortement convexe ? Somme d'une fonction dierentiable et d'une fonction convexe. f

1dierentiable + gradient Lipschitz)

f

1(x)f1(xk) +hrf1(xk);xxki+L2

jjxxkjj22f

1(x)+f2(x)f1(xk) +hrf1(xk);xxki+L2

jjxxkjj22+f2(x) |{z} k(x)

En prenantxk+1= argmin

x2Rnk(x) on assure : f(xk+1)f(xk) | ou|est une \generalisation" dekrf(xk)k22L. Somme d'une fonction dierentiable et d'une fonction convexe. f

1dierentiable + gradient Lipschitz)

f

1(x)f1(xk) +hrf1(xk);xxki+L2

jjxxkjj22f

1(x)+f2(x)f1(xk) +hrf1(xk);xxki+L2

jjxxkjj22+f2(x) |{z} k(x)

En prenantxk+1= argmin

x2Rnk(x) on assure : f(xk+1)f(xk) | ou|est une \generalisation" dekrf(xk)k22L. Somme d'une fonction dierentiable et d'une fonction convexe.

L'iteration precedente est equivalente a :

x k+1= argmin x2Rn1L f2(x) +12 x x krf1(xk)L

Soit encore :

x k+1= proxf2=L x krf1(xk))L Somme d'une fonction dierentiable et d'une fonction convexe.Resultat de convergence 1 [Nesterov 07] Sous les hypotheses precedentes, pourdadmissible tel que jjdjj 1

Df(xk)(d) O

Lf(x0)minfpk

Resultat de convergence 2

Si de plusf1est convexe alors :

f(xk)f(x)Ljjy0yjj2k (variables duales)Nette amelioration par rapport aux techniques de sous-gradient!

Schema multi-pas [Nesterov 07].

In:Nombre d' iterationsN, point initialx02Rn.

Out:xNune estimee dex.

Init:Posert1= 1,y1=x0.

Pourkallant de 0 aN:

Poserxk= proxf2=L

y krf1(yk)L

Calculertk+1=1+p1+4t2k2

Poseryk+1=xk+tk1t

k+1 (xkxk1).

Schema multi-pas [Nesterov 07].

Resultat de convergence

L'algorithme assure que :

f(xk)f(x)Ljjx0xjj2k 2 C'est un taux de convergence optimal.Resultats pratiques... Pour les transformees \simples" de l'analyse harmonique,

30 iterations menent a une precision susante pour le

systeme visuel. Vers du temps reel ? De l'ordre de la seconde pour une image 10001000 (architecture GPU).

Minimisation d'une classe de fonctions fortement

convexes. Le cadre precedent ne permet pas de traiter des problemes de type : min x2Rnkrxk1+12 kxx0k2 en eet : x7! krxk1est trop complexe pour en calculer la resolvante. x7! krxk1n'est pas dierentiable...

On va voir que

la dualit e p ermetde con tournerces probl emes.

Une classe de fonctions fortement convexes.

min x2RnP(x) :=f1(Ax) +f2(x) f1:Rn!R[ f+1gconvexe, s.c.i.

A:Rn!Rmlineaire (ou ane).

f2:Rn!R[ f+1g-fortement convexe. Ari(dom(f2))\ri(dom(f1))6=;.Exemple typique (TV + termel2) : min x2Rnjjrxjj1+2 jjxx0jj22

Complements d'analyse convexe...

Denition (polaire)

Soitf:Rn!R[ f+1gune fonction convexe s.c.i. La

transformee de Fenchel ou polaire defest denie par f (x0) = sup x2Rnhx0;xi f(x)

Proprietes

fest convexe s.c.i.

On af=fce qui implique que :

f(x0) =f(x0) = sup x2Rnhx0;xi f(x)

On peut redenirfa partir de sa polaire.

Complements d'analyse convexe...

Dualite de Fenchel-Rockafellar :On a donc :

min x2Rnf1(Ax) +f2(x) |{z}

Probleme primal

= min x2Rnsup y2RmhAx;yi f1(y) +f2(x) = min x2Rnsup y2Rmhx;Ayi f1(y) +f2(x) = sup y2Rminfx2Rnhx;Ayi f1(y) +f2(x) (Voir schema) = sup y2Rmf1(y)sup x2Rnhx;Ayi f2(x) = sup y2Rmf1(y)f2(Ax) =infy2Rmf1(y) +f2(Ax) |{z}

Probleme dual

Complements d'analyse convexe...

Figure:Une fonction convexe-concave et son point-selle.

Complements d'analyse convexe...

Relations primales-duales.

Soit (x;y) un point selle, on a :

02Ax@f1(y)

02 Ay@f2(x)

Soit encore :

y

2(@f1)(Ax)

x

2(@f2)(Ay)

Car : x

22@f(x1),x12@f(x2)

Les solutions du probleme dual apportent des informations sur le primal et vice-versa !

Complements d'analyse convexe...

Theoreme (dernieres pages du livre de J.B.

Hiriart-Urruty et C. Lemarechal):Soit

f:Rn!R[ f+1gune fonction convexe, s.c.i. Les deux propositions suivantes sont equivalentes :

1.fest dierentiable etrfestL-Lipschitz.

2.fest fortement convexe de module1L

Minimisation d'une classe de fonctions fortement

convexes.Conclusion

Sif2est-fortement convexe, on a :

min x2Rnf1(Ax) +f2(x) =miny2Rmf2(Ay) +f1(y) f2est dierentiablea vecun gradien t1 -Lipschitz. f1est convexe. Pour touty2Y,x=rf2(Ay).Un algorithme naturel : gradient proximal y k+1= proxf1=Lyk+Arf2(Ayk) x k+1=rf2(Ayk+1):

Minimisation d'une classe de fonctions fortement

convexes.Analyse du taux de convergence.

On a :

min x2Rnf1(Ax) +f2(x)|{z}

P(x)=miny2Rmf2(Ay) +f1(y)|{z}

D(y) et pour touty2Y,x=rf2(Ay).Proposition [Peyre, Fadili, Weiss] (sur demande) Soit x(y) =rf2(Ay) Alors kx(y)xk22 jD(y)D(y)j Un taux de CV sur le dual)taux de CV sur le primal.

Minimisation d'une classe de fonctions fortement

convexes.Resultat de convergence (gradient proximal dual)

Sous les hypotheses precedentes :

jjxkxjj2jjjAjjj2jjy0yjj2k En notant (xk;yk) =P(xk)D(yk), le saut de dualite on a : jjxkxjj2(xk;yk):La distance au minimiseur peut ^etre evaluee a chaque iteration !

Minimisation d'une classe de fonctions fortement

convexes..Exemple : le probleme TV-L2 min x2Rnjjrxjj1+2 jjxx0jj2=miny2R2n;jjyjj11jjryx0jj2

Similaire a :

[A. Chambolle 04] An algorithm for TV minimization.

Le schema assure que :

jjxkxjj28nk

La complexite augmente lineairement avec la

dimension !

Minimisation d'une classe de fonctions fortement

convexes.Taux de convergence Les schemas acceleres de Nesterov appliques au dual assurent : jjxkxjj22jjjAjjj2jjy0yjj2 2k2 et sidom(f1) =Rm: f(xk)f(x)2jjy0yjj2jjjAjjj2k 2: (Sinon, rien n'assure que les iterees primales soient admissibles...)

Quelques applications pratiques des algorithmes

precedents

On considere le probleme :

min x2Rnkrxk1+2 kxx0k22 Debruitage dex0sous l'hypothese de bruit gaussien additif. [FILM, AUTRES EXEMPLES]

Quelques applications pratiques des algorithmes

precedents

Debruitage d'images SPIM.

Travail en collaboration avec Jer^ome Fehrenbach (Institut de

Mathematiques de Toulouse) et le canceropole.

L'existence de schemas ecaces doit guider la modelisation.

Motivation et presentation de la methode

Images SPIM:

les images son taect eespar des raies som bres (ou parfois brillantes).

Motivation et presentation de la methode

On se propose d'

eliminerces rai es , pour obtenir des images plus lisibles ãL'origine physique de ces raies est peu claire. On va donc adopter une methode de traitemen td 'image sans modelisation de la physique.

Motivation et presentation de la methode

Notre algorithme s'inspire des methodes recentes de decomposition d'images en :

ãu0=u+v(structure + texture);.

ãu0=u+v+b(structure + texture + bruit):

[Meyer 01, Starck Donoho 05, Aujol 05, Fadili 10, ]=+

Motivation et presentation de la methode

ãDiculte: grand nom bred'inconn ues(plus de 500 images

100010001-50Go.).

ãPour s'en sortir: form ulerle probl emecomme un probl eme de minimisation con vexe . Cela permet le developpement d' algorithmes ecaces , qui convergent vers un minim um global

Modelisation du bruit

On decrit le bruit comme etant compose de raies + bruit blanc (deux composantes).

La raie est un \

motif elementaire " qui se retrouve a plusieurs endroits dans l'image. Pour d eplacer et r epliquer un motif dans l'espace on utilise un pro duitde con volution (x) =X y(y) (xy):

Hypothese sous-jacente :

le bru itest stationnaire

Modelisation du bruit

De facon plus generale, le bruit est modelise comme le produit de convolution entre un bruit blanc (pas forc ementgaussien) et un noyau donne.Figure:Exemples de bruits stationaires obtenus par convolution d'un noyau avec un bruit blanc.

Modelisation du bruit

Autres exemples de produit de convolution:=

Le produit de convolution s'exprime facilement en Fourier.

Modelisation du bruit

Modelisation d'une raie : unltre de Gabor(=une gaussienne anisotrope modulee par une sinusode).

Modelisation du bruit

Le bruit est decrit par plusieurs types de raies (longueur/largeur dierentes). Il s'ecrit b=mX i=1 i i: La loi de probabilite de chaqueidepend de la modelisation du bruit.

D'une maniere generale on ecrit:

P(i) = exp(i(i)):

L'hypothese que le bruit est blanc permet de choisir des fonctionsi separables (ce qui simplie l'analyse n umerique et est justi epar l'hypothese de stationnar ite

Modele de restauration

On utilise le principe du

Maxim umA P osteriori

. On cherche l' imageula plus probable etantdonn ees:

1)la forme et l'orientation dec haqueraie i

(en pratique : une/deux raie + un bruit Gaussien),

2)une description dela probabilit ede c haquebruit i,

3)une description d'unmo deled'image a p riori.

La recherche de l'imageula plus probable peut s'exprimer comme un probl emede minimisation con vexe

Modele de restauration

Principe du Maximum A Posteriori (MAP) :

Connaissantu02Rn, trouver Argmax

2RmnP(ju0):

Loi de Bayes :P(ju0) =P(u0j)P()P(u0);

donc sous l'hypothese d'ind ependancede uet de:

Argmax

2RmnP(ju0) = Argmin

2Rmn(log(P(u))log(P()))

Soit encore :

Argmin

2Rmn0 r u 0mX i=1 ii! 1;+mX i=1 i(i)1 A

Resultats (verication du principe)

Argmin

2R3n r u3X i=1 i? i!

1+k1k1+k2k22+k3k1

Noyaux : \raie verticale (l1)", \Dirac (l2)", "poisson (l1)".=

Resultats (verication du principe)

Resultats (verication du principe)

Figure:Le ltre naturel : un ltre en forme de gradient topologique... Mais mauvaise resolution.Figure:Super-resolution... Le cerveau de l'ombre.

Resultats (verication du principe)

Figure:Le ltre naturel : un ltre en forme de gradient topologique... Mais mauvaise resolution.Figure:Super-resolution... Le cerveau de l'ombre.

Resultats (verication du principe)

Figure:Le ltre naturel : un ltre en forme de gradient topologique... Mais mauvaise resolution.Figure:Super-resolution... Le cerveau de l'ombre.

Resultats

Resultats

Figure:de gauche a droite : originale; 2 ltres (raie de Gabor + Dirac pour le bruit gaussien) ; methode TV-L2

Resultats

Sur une autre image :=+

Remarque sur le choix des fonctionsi

Le modele de decomposition est le suivant :

Argmin

quotesdbs_dbs29.pdfusesText_35

[PDF] Exercices de gestion des ressources humaines-2 - Numilog

[PDF] Mecanique quantique Cours et exercices corriges - Numilog

[PDF] Correction de l 'examen de probabilités et statistiques - Julian Tugaut

[PDF] Examen final corrigé (janvier 2013)

[PDF] Série d 'exercices Finance Internationale - ResearchGate

[PDF] Cours de gestion financière (M1) Exercices corrigés Le cas

[PDF] Exercice n°1

[PDF] Réseaux de Petri

[PDF] Examen professionnel Informatique, système d 'information

[PDF] Graphes exercices et correction

[PDF] Correction examen théorie des jeux 2009-2010 - Ceremade

[PDF] Exercice n° HU 0601 - Corrigé

[PDF] 10

[PDF] 14

[PDF] Examen d analyse personnel technique (ANT) - carrieres gouv