Groupe n˚2 - Equipe n˚4

Projet n˚2
COMPRESSION D’IMAGE A TRAVERS Responsable : xluo
LA FACTORISATION ”SVD” Secrétaire : imhamdi
Codeurs : achasseloupdechatill,wandujar,ygarbage

Résumé : Ce projet consiste à mettre en place un algorithme permettant de faire la compression d’images.
Pour cela on utilisera des techniques matricielles basées sur la factorisation SVD. La première 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éthodes
directe et itérative. La dernière partie applique alors ces transformations à 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ée à un vecteur non nul N de taille n peut s’écrire sous la forme
suivante où id est la matrice identité de dimension n :
N tN
H = Id − 2 × (1)
||N ||2

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ême taille et norme
que X. L’expression mathématique du vecteur N en fonction de X et Y est la suivante :
X −Y
N= (2)
||X − Y ||

Nous avons donc mis en place un algorithme que nous avons appelé obtenirVecteurSymetrie(X,Y)
qui permet à l’aide de l’équation (2) d’obtenir le vecteur de symétrie N définissant la matrice de
HouseHolder. Ensuite, pour obtenir la matrice de HouseHolder à partir des vecteurs X et Y nous
avons défini la fonction obtenirHouseHolderXY(X, Y) qui à l’aide de l’algorithme précédent permet
de trouver cette matrice en utilisant l’équation (1). L’algorithme obtenirHouseHolderN(N) fait la
même chose mais il prend comme paramètre le vecteur de symétrie N .

     
3 0 0.64 −0.48 0.6
Ainsi pour envoyer 4 sur 0 la matrice de HouseHolder est −0.48 0.36 0.8
0 5 0.6 0.8 0

1
1.2 Optimisation du produit matrice-matrice
La multiplication d’une matrice de HouseHolder H par un vecteur U revient simplement à trouver
l’image de ce vecteur par la matrice de HouseHolder en utilisant l’expression initiale de H c’est à dire :

N tN N tN
 
H × U = Id − 2 × U = U − 2 × U (3)
||N ||2 ||N ||2

Pour cela nous avons créé une fonction appliquerHouseHolderVecteur prenant comme paramètres
les deux vecteurs X et Y définissant notre matrice de HouseHolder et le vecteur U . Cet algorithme
calcule le produit matriciel à l’aide de la fonction houseHolderInterne traduisant l’équation (3).
Nous souhaitons optimiser le produit matriciel entre une matrice de HouseHolder et une autre ma-
trice de même taille. Pour éviter de faire un produit matrice-matrice nous avons utilisé la fonction
précédente 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ée de N
L ← nombre de lignes de N
C ← nombre de colonnes de M
Res ← matrice L × C

pour i de 0 à L faire
iième ligne de Res ← houseHolderInterne(N, tN, iième colonne de M)
retourner Res

La complexité de cette fonction est de l’ordre de n produits matrice-vecteur (où n est la taille de
la matrice). D’autre part, d’après l’algorithme appliquerHouseHolderVecteur un produit matrice-
vecteur est équivalent à deux produits vecteur-vecteur donc de complexité 2 × n. Finalement, la
complexité de notre algorithme est de l’ordre de 2 × n × n équivalent à Θ(n2 ). Alors que le produit
matrice-matrice usuel admet une complexité de l’ordre de n × n × Θ(vecteur − vecteur) c’est-à-dire
une complexité équivalente à Θ(n3 ).

2 Mise sous forme bidiagonale

Dans cette partie nous expliquons les algorithmes permettant de mettre une matrice carrée sous
forme bidiagonale. Pour ce faire, nous devons extraire des vecteurs colonnes et des vecteurs lignes
d’une matrice donnée. Nous avons également besoin de générer des vecteurs à une seule composante
non nulle.

2.1 Extraction de vecteurs colonnes et lignes

Nous avons codé deux fonctions extraireVecteurColonne(M, i) (respectivement la fonction
extraireVecteurLigne(M, i)) permettant d’extraire le vecteur colonne de la ième colonne (respec-
tivement le vecteur ligne de la iè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 fonc-
tions 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 1ère ligne et la 1ère colonne de la matrice A. En effet,
nous avons respecté les contraintes de l’algorithme de bidiagonalisation suivantes :
– le vecteur extrait de la ième colonne a ses i premières composantes nulles,
– le vecteur extrait de la ième ligne a ses i+1 premières composantes nulles.

2.2 Génération de vecteurs de même norme à une seule composante non

nulle
Nous avons également mis en place deux fonctions genererVecteurColonne(M, i) (respective-
ment genererVecteurLigne(M, i)) permettant de générer le vecteur de même norme que le vecteur
colonne extrait de la ième colonne (respectivement le vecteur ligne extrait de la ième ligne) et ayant
seulement sa ième (respectivement (i + 1)ème ) composante non nulle.
Ainsi, en considérant la même matrice A que précédemment, 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ée

Ces fonctions auxiliaires étant disponibles, nous avons implémenté l’algorithme de mise sous forme
bidiagonale suivant :

fonction miseSousFormeBidiagonale(A : matrice carrée)

n ← taille(A)
Qlef t , Qright ← matriceIdentite(n)
BD ← A
pour i de 0 à 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ée, ainsi que les changements de base à gauche et à droite.
À chaque tour de boucle de l’algorithme, nous pouvons vérifier que Qlef t · BD · Qright = 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 écrit une fonction SVD decompositionQR qui prend comme
paramètres une matrice bidiagonale BD et un entier représentant le nombre d’itérations à effectuer
et qui renvoie un triplet de matrices U , S et V ayant des propriétés spécifiques. Cet algorithme utilise
la fonction [Link] 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égligeant les termes < 10−5 .
S=  0 0 1.73675226 0 
0 0 0 0.64329448
Pour montrer cette propriété nous avons mis en place des fonctions permettant de tester la conver-
gence de S. La fonction dessinerConvergenceDiagS permet de dessiner cette convergence. Elle prend
comme paramètres la taille de la matrice BD, la précision et le nombre d’itérations. L’algorithme est
le suivant :

fonction dessinerConvergenceDiagS (n, step : entiers, precision : réel)

BD ← matrice bidiagonale aléatoire de taille n
pour i de 1 à step faire
res ← SVD decompositionQR (BD, i)
S ← 2éme élément de la liste res
y[i] ← nombre d’éléments extradiagonaux de S
Tracer le graphique y = f (step)

Pour une matrice bidiagonale de taille n et pour un nombre d’itération égal à 100 nous avons obtenu
la Figure 1 (voir page suivante).

D’autre part, nous avons créé une fonction verifierInvariantEgalité USV BD pour s’assurer que
l’invariant U × S × V = BD est toujours vérifié (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érations

3.2 Optimisation de la décomposition 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ètre une matrice bidi-
agonale. Donc, S est une matrice bidigonale dès le départ. Or, d’après 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érieure
⇒ t S est une matrice bidiagonale inférieure.
⇒ 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érieure (car R1 est une matrice triangulaire supérieure - d’après
la décomposition QR - et tridiagonale).
On peut suivre le même raisonnement pour R2 pour montrer l’invariant.

Nous allons maintenant calculer la complexité de l’algorithme de factorisation SVD sachant que
nous n’avons pas réussi à implémenter l’algorithme optimisé de la factorisation QR. Notons Θ(QR) la
complexité de la factorisation QR. La complexité de l’algorithme de décomposition SVD est de l’ordre
de step × Θ(QR) or d’après le cours la complexité de l’algorithme de décomposition QR est de l’ordre
de step × Θ(n2 ) et peut même aller jusqu’à step × Θ(n) si la matrice de départ est bidigonale. Donc
la complexité de la fonction de factorisation SVD optimisée est de l’ordre de Θ(step2 × n).

5
3.3 Modification des matrices U et S
La décomposition SVD demande à ce que les éléments de la matrice S soient positifs et ordonnés
de manière décroissante. Nous avons donc créé une fonction modifierUS qui prend comme paramètre
les matrices U et S et les modifie de sorte que S vérifie la propriété précédente et que l’invariant
(U × S × V ) = BD soit conservé.

fonction modifierUS (U, S : matrice)

l ← nombre de lignes de S
C ← nombre de colonnes de S
diagS ← diagonale de S
triS ← les éléments de la diagonale de S positifs et triés
diagS ← diagS
triS
pour i de 0 à min(l, C) faire
S[i,i] ← triS[i]
iième colonne de U ← iième colonne de U × diagS[i]
retourner (U,S)

4 Application à la compression d’image

Dans cette partie nous appliquons la transformation SVD à l’algorithme de compression d’image.
Nous avons commencé par mettre en place un algorithme CompressionRangK permettant d’obtenir
une compression au rang k d’une matrice carrée donnée. En effet, cet algorithme modifie le rang de
la matrice S en la copiant jusqu’à la colonne k et en mettant à 0 le reste des valeurs.

4.1 Gain pour la compression au rang k

Soit A la matrice carré de dimension n représentant l’image d’origine. Notons A0 la matrice de
l’image compressée. En décomposant les matrices U , S et V par blocs, on obtient le calcul suivant :
     
0 U1 U2 S1 0 V1 V2 U1 S1 V1 U1 S1 V2
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ésultat final. Pour
stocker la matrice A0 , il suffit donc de mémoriser 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 stricte-
ment négatif, la transformation est de taille supérieure à l’image initiale. En revanche, pour que la
n2
compression soit efficace, k ne doit pas dépasser 2n+1 .

4.2 Efficacité de la compression

Afin de pouvoir déterminer l’efficacité, il a fallu utiliser un exemple. Nous avons choisi l’exemple
de l’image de la NASA représentant la terre (500pxpar500px). Pour pouvoir l’exploiter, il nous a fallu
réaliser son acquisition, son traitement et sa recomposition.

6
4.2.1 Acquisition, traitement et recomposition de l’image
L’image est encodée sous forme d’une matrice de pixels, eux même représentés sur la base du
système RGB par des matrices de trois éléments. Une fois l’acquisition effectuée, il faut décomposer
cette matrice en trois matrices afin de pouvoir travailler sur chaque couleur des pixels. Pour ce faire,
nous avons utilisé la fonction suivante :
fonction decomposerMatrice(M : Matrice)
n ← taille(S)
M1, M2, m3 ← CreerMatriceCarre(n)

pour i de 0 à k faire
pour j de k à 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écomposition 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ée sur le même principe que l’algorithme de décomposition, nous
ne détaillerons 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 à gauche l’image d’origine.

image initiale compression au rang 50 compression au rang 100

Figure 2 – compression d’une image au rang 50 et 100

7
4.2.3 Efficacité
L’efficacité de la compression se détermine à partir de la distance entre la matrice de départ et
la matrice de l’image compressée. Pour cela, on utilise la norme fournie par la bibliothèque numpy
appliquée sur la différence entre les deux matrices.

Efficacité(Mat dep, Mat res) = ||Mat dep - Mat res ||

Ainsi, en itérant cette méthode pour des compressions allant de k = 0 à k = 500 on obtient la
courbe de l’efficacité suivante :

Figure 3 – Courbe de l’efficacité de la compression en fonction du rang k

5 Conclusion
Ce projet nous a paru particulièrement intéressant pour ce qui est de la manipulation des matrices
de HouseHolder et de l’implémentation 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.