Transcript 1 Transformations de Householder

Projet n˚2
COMPRESSION D’IMAGE A TRAVERS
LA FACTORISATION ”SVD”
Groupe n˚2 - Equipe n˚4
Responsable : xluo
Secr´etaire : imhamdi
Codeurs : achasseloupdechatill,wandujar,ygarbage
R´
esum´
e : Ce projet consiste `
a mettre en place un algorithme permettant de faire la compression d’images.
Pour cela on utilisera des techniques matricielles bas´ees sur la factorisation SVD. La premi`ere partie de ce
projet met en place les fonctions de manipulation des matrices de HouseHolder qu’on utilisera dans les deux
parties suivantes afin de transformer une matrice A en une matrice bidiagonale BD avec les deux m´ethodes
directe et it´erative. La derni`ere partie applique alors ces transformations `
a l’algorithme de compression
d’image par transformation SVD.
1
Transformations de Householder
Dans cette partie, nous allons expliquer les algorithmes utilitaires mis en place pour manipuler les
matrices de Householder.
1.1
La matrice de HouseHolder
Une matrice de HouseHolder associ´ee `a un vecteur non nul N de taille n peut s’´ecrire sous la forme
suivante o`
u id est la matrice identit´e de dimension n :
H = Id − 2 ×
N tN
||N ||2
(1)
Nous souhaitons dans un premier temps trouver l’expression du vecteur N de taille n pour que la
matrice de HouselHolder H envoie le vecteur X (de taille n) sur le vecteur Y de mˆeme taille et norme
que X. L’expression math´ematique du vecteur N en fonction de X et Y est la suivante :
N=
X −Y
||X − Y ||
(2)
Nous avons donc mis en place un algorithme que nous avons appel´e obtenirVecteurSymetrie(X,Y)
qui permet `
a l’aide de l’´equation (2) d’obtenir le vecteur de sym´etrie N d´efinissant la matrice de
HouseHolder. Ensuite, pour obtenir la matrice de HouseHolder `a partir des vecteurs X et Y nous
avons d´efini la fonction obtenirHouseHolderXY(X, Y) qui `a l’aide de l’algorithme pr´ec´edent permet
de trouver cette matrice en utilisant l’´equation (1). L’algorithme obtenirHouseHolderN(N) fait la
mˆeme chose mais il prend comme param`etre le vecteur de sym´etrie N .
 
 

3
0
0.64 −0.48
Ainsi pour envoyer 4 sur 0 la matrice de HouseHolder est −0.48 0.36
0
5
0.6
0.8
1

0.6
0.8
0
1.2
Optimisation du produit matrice-matrice
La multiplication d’une matrice de HouseHolder H par un vecteur U revient simplement `a trouver
l’image de ce vecteur par la matrice de HouseHolder en utilisant l’expression initiale de H c’est `a dire :
N tN
N tN
U
=
U
−
2
×
U
(3)
H × U = Id − 2 ×
||N ||2
||N ||2
Pour cela nous avons cr´e´e une fonction appliquerHouseHolderVecteur prenant comme param`etres
les deux vecteurs X et Y d´efinissant notre matrice de HouseHolder et le vecteur U . Cet algorithme
calcule le produit matriciel `
a l’aide de la fonction houseHolderInterne traduisant l’´equation (3).
Nous souhaitons optimiser le produit matriciel entre une matrice de HouseHolder et une autre matrice de mˆeme taille. Pour ´eviter de faire un produit matrice-matrice nous avons utilis´e la fonction
pr´ec´edente qui fait intervenir uniquement des produits matrice-vecteur. L’algorithme mis en place est
le suivant :
fonction appliquerHouseHolderMatrice(X,Y : vecteurs, M : matrice)
N ← obtenirVecteurSymetrie(X, Y)
tN ← transpos´
ee de N
L ← nombre de lignes de N
C ← nombre de colonnes de M
Res ← matrice L × C
pour i de 0 `
a L faire
ii`eme ligne de Res ← houseHolderInterne(N, tN, ii`eme colonne de M)
retourner Res
La complexit´e de cette fonction est de l’ordre de n produits matrice-vecteur (o`
u n est la taille de
la matrice). D’autre part, d’apr`es l’algorithme appliquerHouseHolderVecteur un produit matricevecteur est ´equivalent `
a deux produits vecteur-vecteur donc de complexit´e 2 × n. Finalement, la
complexit´e de notre algorithme est de l’ordre de 2 × n × n ´equivalent `a Θ(n2 ). Alors que le produit
matrice-matrice usuel admet une complexit´e de l’ordre de n × n × Θ(vecteur − vecteur) c’est-`a-dire
une complexit´e ´equivalente `
a Θ(n3 ).
2
Mise sous forme bidiagonale
Dans cette partie nous expliquons les algorithmes permettant de mettre une matrice carr´ee sous
forme bidiagonale. Pour ce faire, nous devons extraire des vecteurs colonnes et des vecteurs lignes
d’une matrice donn´ee. Nous avons ´egalement besoin de g´en´erer des vecteurs `a une seule composante
non nulle.
2.1
Extraction de vecteurs colonnes et lignes
Nous avons cod´e deux fonctions extraireVecteurColonne(M, i) (respectivement la fonction
extraireVecteurLigne(M, i)) permettant d’extraire le vecteur colonne de la ie`me colonne (respectivement le vecteur ligne de la ie`me ligne) de la matrice M . Cependant, pour pouvoir utiliser les
2
fonctions fournies par la partie 1 permettant d’obtenir des matrices de HouseHolder, les deux fonctions d’extraction produisent des vecteurs colonnes.


11 21 5 12
48 548 16 887

Exemple : si on se donne la matrice A suivante : 
47 88 91 45 , les appels des
17 32 68 17
fonctions extraireVecteurColonne(A, 1) et extraireVecteurLigne(A,
1) produisent les
vecteurs : t 0 58 88 32 et t 0 0 16 887
Les vecteurs ci-dessus ne sont pas exactement la 1e`re ligne et la 1e`re colonne de la matrice A. En effet,
nous avons respect´e les contraintes de l’algorithme de bidiagonalisation suivantes :
– le vecteur extrait de la ie`me colonne a ses i premi`eres composantes nulles,
– le vecteur extrait de la ie`me ligne a ses i+1 premi`eres composantes nulles.
2.2
G´
en´
eration de vecteurs de mˆ
eme norme `
a une seule composante non
nulle
Nous avons ´egalement mis en place deux fonctions genererVecteurColonne(M, i) (respectivement genererVecteurLigne(M, i)) permettant de g´en´erer le vecteur de mˆeme norme que le vecteur
colonne extrait de la ie`me colonne (respectivement le vecteur ligne extrait de la ie`me ligne) et ayant
seulement sa ie`me (respectivement (i + 1)e`me ) composante non nulle.
Ainsi, en consid´erant la mˆeme matrice A que pr´ec´edemment, les appels de genererVecteurColonne(A,
1) et genererVecteurLigne(A, 1) produisent les vecteurs : t 0 556 0 0 et t 0 0 887.1 0
2.3
Mise sous forme bidiagonale d’une matrice carr´
ee
Ces fonctions auxiliaires ´etant disponibles, nous avons impl´ement´e l’algorithme de mise sous forme
bidiagonale suivant :
fonction miseSousFormeBidiagonale(A : matrice carr´
ee)
n ← taille(A)
Qlef t , Qright ← matriceIdentite(n)
BD ← A
pour i de 0 `
a n − 1 faire
X ← extraireVecteurColonne(BD, i)
singleX ← genererVecteurColonne(BD, i)
Q1 ← obtenirHouseHolderXY(X, singleX)
Qlef t ← Qlef t · Q1
BD ← Q1 · BD
si (i < (n − 2)) alors
X ← extraireVecteurLigne(BD, i)
singleX ← genererVecteurLigne(BD, i)
Q2 ← obtenirHouseHolderXY(X, singleX) ;
Qright ← Q2 · Qright
BD ← BD · Q2
retourner (Qlef t , BD, Qright )
3
L’algorithme renvoie la matrice transform´ee, ainsi que les changements de base `a gauche et `a droite.
` chaque tour de boucle de l’algorithme, nous pouvons v´erifier que Qlef t · BD · Qright = A.
A
3
Transformations QR
Dans cette partie nous allons mettre sous forme SVD une matrice diagonale en appliquant un
certain nombre de fois la transformation QR.
3.1
Factorisation SVD
Dans un premier temps nous avons ´ecrit une fonction SVD decompositionQR qui prend comme
param`etres une matrice bidiagonale BD et un entier repr´esentant le nombre d’it´erations `a effectuer
et qui renvoie un triplet de matrices U , S et V ayant des propri´et´es sp´ecifiques. Cet algorithme utilise
la fonction numpy.linalg.qr pour la transformation QR.
Dans cet algorithme la matrice S converge vers une matrice diagonale :


1 4 0 0
0 4 1 0

Exemple : L’application de l’algorithme sur la matrice 
0 0 3 4 donne la matrice
0 0 0 3


5.85636553
0
0
0


0
5.50207658
0
0
 en n´
egligeant les termes < 10−5 .
S=


0
0
1.73675226
0
0
0
0
0.64329448
Pour montrer cette propri´et´e nous avons mis en place des fonctions permettant de tester la convergence de S. La fonction dessinerConvergenceDiagS permet de dessiner cette convergence. Elle prend
comme param`etres la taille de la matrice BD, la pr´ecision et le nombre d’it´erations. L’algorithme est
le suivant :
fonction dessinerConvergenceDiagS (n, step : entiers, precision : r´
eel)
BD ← matrice bidiagonale al´
eatoire de taille n
pour i de 1 `
a step faire
res ← SVD decompositionQR (BD, i)
S ← 2e´me ´
el´
ement de la liste res
y[i] ← nombre d’´
el´
ements extradiagonaux de S
Tracer le graphique y = f (step)
Pour une matrice bidiagonale de taille n et pour un nombre d’it´eration ´egal `a 100 nous avons obtenu
la Figure 1 (voir page suivante).
D’autre part, nous avons cr´e´e une fonction verifierInvariantEgalit´
e USV BD pour s’assurer que
l’invariant U × S × V = BD est toujours v´erifi´e (ce qui est le cas voir tests de la partie 3 ).
4
Figure 1 – Convergence de S vers une matrice bidiagonale en fonction du nombre d’it´erations
3.2
Optimisation de la d´
ecomposition QR
Nous allons maintenant montrer l’invariant S,R1 et R2 sont des matrices bidiagonales.
L’algorithme de factorisation SVD SVD decompositionQR prend comme param`etre une matrice bidiagonale. Donc, S est une matrice bidigonale d`es le d´epart. Or, d’apr`es le cours la factorisation QR
conserve la forme d’une matrice tridiagonale. Pour montrer l’invariant nous pouvons alors suivre le
raisonnement suivant :
S est une matrice bidiagonale sup´erieure
⇒ t S est une matrice bidiagonale inf´erieure.
⇒ t S est une matrice tridiagonale.
⇒ R1 est une matrice tridiagonale (car la factorisation QR conserve la forme tridiagonale).
⇒ R1 est une matrice bidiagonale sup´erieure (car R1 est une matrice triangulaire sup´erieure - d’apr`es
la d´ecomposition QR - et tridiagonale).
On peut suivre le mˆeme raisonnement pour R2 pour montrer l’invariant.
Nous allons maintenant calculer la complexit´e de l’algorithme de factorisation SVD sachant que
nous n’avons pas r´eussi `
a impl´ementer l’algorithme optimis´e de la factorisation QR. Notons Θ(QR) la
complexit´e de la factorisation QR. La complexit´e de l’algorithme de d´ecomposition SVD est de l’ordre
de step × Θ(QR) or d’apr`es le cours la complexit´e de l’algorithme de d´ecomposition QR est de l’ordre
de step × Θ(n2 ) et peut mˆeme aller jusqu’`a step × Θ(n) si la matrice de d´epart est bidigonale. Donc
la complexit´e de la fonction de factorisation SVD optimis´ee est de l’ordre de Θ(step2 × n).
5
3.3
Modification des matrices U et S
La d´ecomposition SVD demande `
a ce que les ´el´ements de la matrice S soient positifs et ordonn´es
de mani`ere d´ecroissante. Nous avons donc cr´e´e une fonction modifierUS qui prend comme param`etre
les matrices U et S et les modifie de sorte que S v´erifie la propri´et´e pr´ec´edente et que l’invariant
(U × S × V ) = BD soit conserv´e.
fonction modifierUS (U, S : matrice)
l ← nombre de lignes de S
C ← nombre de colonnes de S
diagS ← diagonale de S
triS ← les ´
el´
ements de la diagonale de S positifs et tri´
es
diagS ← diagS
triS
pour i de 0 `
a min(l, C) faire
S[i,i] ← triS[i]
ii`eme colonne de U ← ii`eme colonne de U × diagS[i]
retourner (U,S)
4
Application `
a la compression d’image
Dans cette partie nous appliquons la transformation SVD `a l’algorithme de compression d’image.
Nous avons commenc´e par mettre en place un algorithme CompressionRangK permettant d’obtenir
une compression au rang k d’une matrice carr´ee donn´ee. En effet, cet algorithme modifie le rang de
la matrice S en la copiant jusqu’`
a la colonne k et en mettant `a 0 le reste des valeurs.
4.1
Gain pour la compression au rang k
Soit A la matrice carr´e de dimension n repr´esentant l’image d’origine. Notons A0 la matrice de
l’image compress´ee. En d´ecomposant les matrices U , S et V par blocs, on obtient le calcul suivant :
U1 U2
S1 0
V1 V2
U1 S1 V1 U1 S1 V2
0
A =
=
(4)
U3 U4
0 0
V3 V4
U3 S1 V1 U3 S1 V2
On remarque alors que les blocs U2 , U4 ,V3 et V4 n’interviennent pas dans le r´esultat final. Pour
stocker la matrice A0 , il suffit donc de m´emoriser les valeurs des blocs S1 , U1 , U3 , V1 , et V2 , soit
k + 2k(n − k) + 2k 2 = 2kn + k valeurs.
Le gain de place pour une compression au rang k vaut donc n2 − 2kn − k. Si ce gain est strictement n´egatif, la transformation est de taille sup´erieure `a l’image initiale. En revanche, pour que la
n2
.
compression soit efficace, k ne doit pas d´epasser 2n+1
4.2
Efficacit´
e de la compression
Afin de pouvoir d´eterminer l’efficacit´e, il a fallu utiliser un exemple. Nous avons choisi l’exemple
de l’image de la NASA repr´esentant la terre (500pxpar500px). Pour pouvoir l’exploiter, il nous a fallu
r´ealiser son acquisition, son traitement et sa recomposition.
6
4.2.1
Acquisition, traitement et recomposition de l’image
L’image est encod´ee sous forme d’une matrice de pixels, eux mˆeme repr´esent´es sur la base du
syst`eme RGB par des matrices de trois ´el´ements. Une fois l’acquisition effectu´ee, il faut d´ecomposer
cette matrice en trois matrices afin de pouvoir travailler sur chaque couleur des pixels. Pour ce faire,
nous avons utilis´e la fonction suivante :
fonction decomposerMatrice(M : Matrice)
n ← taille(S)
M1, M2, m3 ← CreerMatriceCarre(n)
pour i de 0 `
a k faire
pour j de k `
a n faire
M1[i][j] = M[i][j][0]
M2[i][j] = M[i][j][1]
M3[i][j] = M[i][j][2]
retourner M1, M2, M3
Une fois l’image acquise, on applique l’algorithme de compression en utilisant la d´ecomposition SVD,
puis en appliquant l’algorithme de compression au rang k pour chacune des trois matrices issues de
la SVD.
La recomposition de l’image est bas´ee sur le mˆeme principe que l’algorithme de d´ecomposition, nous
ne d´etaillerons pas son contenu.
4.2.2
Exemples
Voici ci-dessous deux exemples d’application de nos algorithmes de compression au rang 50 et au
rang 100 avec `
a gauche l’image d’origine.
image initiale
compression au rang 50
Figure 2 – compression d’une image au rang 50 et 100
7
compression au rang 100
4.2.3
Efficacit´
e
L’efficacit´e de la compression se d´etermine `a partir de la distance entre la matrice de d´epart et
la matrice de l’image compress´ee. Pour cela, on utilise la norme fournie par la biblioth`eque numpy
appliqu´ee sur la diff´erence entre les deux matrices.
Efficacit´
e(Mat dep, Mat res) = ||Mat dep - Mat res ||
Ainsi, en it´erant cette m´ethode pour des compressions allant de k = 0 `a k = 500 on obtient la
courbe de l’efficacit´e suivante :
Figure 3 – Courbe de l’efficacit´e de la compression en fonction du rang k
5
Conclusion
Ce projet nous a paru particuli`erement int´eressant pour ce qui est de la manipulation des matrices
de HouseHolder et de l’impl´ementation d’algorithmes de factorisation matricielle (factorisation QR).
Il nous a permis d’apprendre la factorisation SVD et de voir un exemple concret de son application :
la compression d’image.
8