Le premier but de ce projet est d'implémenter astucieusement des algorithmes rapides
de multiplication de matrices. Le premier est dû à Straussen et son amélioration
à Winograd.
Si la multiplication naïve de matrices carrées de taille N requiert
O(N3) opérations (i.e. additions et multiplications), il est possible d'utiliser
seulement O(N2.81) opérations par une technique de division pour régner.
En poursuivant dans cette direction, il serait possible de considérer
des algorithmes plus sophistiqués n'utilisant que O(N2.38) opérations,
mais cela est en dehors de ce projet.
Ensuite, nous aborderons des problèmes connexes utilisant
en sous-routine la multiplication de matrices et
parfois des techniques aléatoires.
Le premier consiste à vérifier qu'un produit de matrice est correct avec seulement O(N2)
opérations. Il s'agit du test probabiliste de Freivalds.
Le deuxième est le calcul de l'inverse d'une matrice
à l'aide de nos algorithmes de multiplication de matrices.
En en déduira un algorithme pour tester si un déterminant est nul.
Une application est
la détection d'un couplage parfait dans un graphe biparti.
Cet application est étroitement reliée au théorème de Schwartz sur le nombre de zéros d'un
polynômes à plusieurs variables dans un corps fini.
Cet algorithme sera comparé à celui de Ford-Fulkerson qui trouve un couplage maximal
en utilisant une technique de maximisation de flot.
Enfin, nous résoudrons le problème de détection de k-cycles dans un graphe non orienté
toujours à l'aide de nos algorithmes de multiplication de matrices.
L'interface n'est pas un objectif dans ce projet.
Toutefois, pour la saisie des matrices ou du graphe (biparti ou non),
le choix sera laissé à l'utilisateur entre un tirage aléatoire, la lecture d'un fichier, ou une saisie manuelle.
Dans le cas d'un fichier, la première ligne du fichier contiendra la taille des matrices,
puis l'énumération des lignes de chaque matrice.
Pour les graphes, en première ligne sera le nombre de noeuds considérés, puis
sur chaque ligne du fichier la liste des successeurs de chaque noeud.
On supposera toujours que les matrices sont de taille N=2n.
Dans votre projet vous pourrez considérer ou expliquer comment considérer
le cas général.
A part pour le calcul de l'inverse d'une matrice, on pourra toujours
se placer sur les entiers afin d'accélérer le temps de calcul et réduire
l'espace mémoire.
Les algorithmes de multiplication rapide sont basés sur une observation
dans le cas des matrices 2× 2, qui se généralise aisément par
induction lorsqu'on découpe une matrice M de taille N en 4 sous-blocs de taille N/2 :
M=
⎛
⎝
M11
M12
M21
M22
⎞
⎠
.
Voici maintenant les décompositions récursivent qu'utilisent respectivement
les algorithmes de Strassen [Str69] et Winograd [Win71].
On ne manquera pas de les justifier.
La matrice C est le produit de deux matrices A et B. Les autres matrices
sont des matrices intermédiaires utilisées dans le calcul.
Pour Strassen :
M1
=(A11 − A22)*(B21 + B22)
M2
=(A11 + A22)*(B11 + B22)
M3
=(A11 − A21)*(B11 + B12)
M4
=(A11 + A12)*B22
M5
=A11*(B12 − B22)
M6
=A22*(B21 − B11)
M7
=(A21 + A22)*B11
C11
=M1 + M2 − M4 + M6
C12
=M4 + M5
C21
=M6 + M7
C22
=M2 − M3 + M5 − M7
Pour Winograd :
S1
=A21 + A22
S2
=S1 − A11
S3
=A11 − A21
S4
=A12 − S2
S5
=B12 − B11
S6
=B22 − S5
S7
=B22 − B12
S8
=S6 − B21
M1
=S2*S6
M2
=A11*B11
M3
=A12*B21
M4
=S3*S7
M5
=S1*S5
M6
=S4*B22
M7
=A22*S8
T1
=M1 + M2
T2
=T1 + M4
C11
=M2 + M3
C12
=T1 + M5 + M6
C21
=T2 − M7
C22
=T2 + M5
Dans un premier temps, la récursion s'arrêtera lorsque les matrices sont de taille 1.
On comparera ces algorithmes avec l'algorithme standard. La vérification proprement
dite du calcul s'effectuera avec la stratégie définie dans le paragraphe suivant.
Cette comparaison mettra en évidence un seuil pour chacun des algorithmes
sur la taille des matrices à partir duquel la multiplication standard est moins rapide.
En dessous, de ce seuil on arêtera la récursion et on calculera de manière standard
le produit.
Freivalds [Fre79] a montré que si deux matrices A et B étaient
différentes, alors AX≠ BX avec probabilité au moins 1/2 lorsque X
est un vecteur aléatoire dont chaque coordonnée est soit 0, soit 1
de manière équiprobable.
Pour vérifier que la matrice C=A× B, il suffit de vérifier
pour un tel X que CX=A(BX), pour suffisamment de vecteurs X.
Vous justifierez que le nombre d'opérations utilisées pour chaque vérification
est en O(N2).
Si C=A× B, alors on a toujours CX=A(BX), si
C=A× B, et qu'on répète T fois cette vérification avec T vecteurs tirés au sort,
alors la probabilité d'avoir CX≠ A(BX) au moins une fois est
au moins 1−1/2T. La probabilité de se tromper est donc de 1/2T.
Dans cette partie il conviendra de se placer sur les réels.
Supposons A inversible, alors A−1=(ATA)−1AT.
Donc il est suffisant de savoir inverser les matrices symétriques définies positives.
Si A est ainsi, alors A−1 se déduit de A par induction sur son découpage en blocs ainsi :
A=
⎛
⎝
B
CT
C
D
⎞
⎠
⇒
A−1=
⎛
⎝
B−1+B−1CTS−1CB−1
−B−1CTS−1
−S−1CB−1
S−1
⎞
⎠
,
avec S=D−CB−1CT.
Une autre manière de procéder à la Strassen est la suivante, où dans ce qui
suit B=A−1 :
P
=(A11)−1
Q
=A21*P
R
=P*A12
S
=A21*R
T
=S−A22
U
=T−1
B12
=R*U
B21
=U*Q
B22
=−U
V
=R*B21
B11
=P−V
On testera chacun des algorithmes implémentés en adaptant le test de Freivalds exposé précédemment.
Puis on adaptera ces algorithmes pour détecter si une matrice est inversible,
i.e. si son déterminant est non nul.
Un graphe equibiparti est la donnée d'un triplet G=(A,B,E),
où A et B sont des ensembles de sommets disjoints et de tailles égales, et E un ensemble
d'arêtes entre A et B. Nous supposerons ici que les arêtes ne sont pas orientées.
Un couplage parfaitC est un sous-ensemble de E qui marie chaque sommet
de A avec un seul sommet de B.
On présente ici l'approche de Mulmuley, Vazirani et Vazirani [MVV87].
On considère des variables xij, pour 1≤ i,j≤ N.
Considérons la matrice M de taille N× N telle que Mij=xij si l'arête (i,j)
est dans E, et Mij=0 sinon.
Le déterminant detM de M est donc un polynôme en (xij) de degré N.
Ce polynôme est formé des monômes de la forme
x1σ(1)x2σ(2)⋯ xNσ(N), où σ est une permutation
de {1,2,…,N} quelconque. Plus précisément ce monôme apparaît dans detM
si et seulement si chacune des arêtes (1,σ(1)),(2,σ(2)),…,(N,σ(N)) est
dans E, qui est alors un couplage parfait du graphe.
Donc detM n'est pas le polynôme nul si et seulement si G admet un couplage parfait.
Le lemme de Schwartz est alors ici utilisé pour vérifier algorithmiquement et aléatoirement
que ce polynôme n'est pas nul.
Soit P un polynôme non nul de degré d à m variables (sur un corps quelconque).
Ce théorème dit que si à chaque variable est assignée à une valeur prise
uniformément au hasard dans un ensemble S de taille s, alors
P ne s'annule pas en ces valeurs avec probabilité au moins d/s.
On comparera les résultats obtenus avec l'algorithme usuel pour calculer
le couplage parfait d'un graphe biparti de Ford-Fulkerson.
Cette méthode utilise la notion de flots pour calculer un couplage de taille
maximale. S'il contient tous les noeuds alors il est parfait.
Pour une description de l'algorithme, on pourra se reporter à la pâle HC 2003-04 :
http://www.enseignement.polytechnique.fr/informatique/IF/04/hc/hc.html.
Si A est la matrice d'adjacence d'un graphe (non orienté),
le calcul de Ak sur l'algèbre booléenne permet de déterminer les
paires de sommets (i,j) reliés par un chemin de taille plus petite que k:
(Ak)i,j=1 ssi un tel chemin existe entre i et j.
Par contre, si A est la matrice d'adjacence d'un graphe orienté acyclique (i.e.
sans cycles), alors (Ak)i,j=1 ssi un chemin existe de longueur
exactement k de i à j.
Etant donné un paramètre fixe k et un graphe non orienté G,
on aimerait déterminer un cycle simple (i.e. sans sous-cycles)
de taille exactement k, et en trouver un s'il en existe,
en utilisant l'approche d'Alon, Yuster et Zwick [AYZ95].
Pour ce, on va orienter aléatoirement G de manière à le rendre acyclique.
Soit le graphe G orienté.
Puis il suffira de vérifier
s'il existe une paire (i,j) connectée par un chemin de longueur (k−1) de i à j dans ,
telle que (i,j) est une arête du graphe initial G non orienté.
L'orientation est effectuée en tirant au sort une permutation σ de V dans {1,2,…, N}, où N
est la taille de V. Alors, une arête (i,j) de G sera une arête de
de i vers j si σ(i)<σ(j), et de j vers i sinon.
On peut montrer que la probabilité de trouver un cycle, s'il en existe, est de 1/k!. Il suffit donc
de repéter T fois l'expérience, en tirant une nouvelle permutation à chaque itération,
pour obtenir une probabilité d'erreur inférieure à (1−1/k!)T.
Le calcul de Ak se fera en O(logk) multiplications de matrices en utilisant
l'algorithme d'exponentiation rapide standard. Même s'il n'est pas possible d'adapter
les algorithmes de multiplication rapide de Strassen et Winograd, il suffira de
faire les multiplications sur les entiers et d'associer chaque coefficient non nul
au booléen 1.
La représentation des matrices est fondamentale dans ce projet.
Une représentation de type array n'est pas adaptée à la
récursivité des méthodes utilisées.
Une représentation récursive en forme d'arbre de rang 4 (chaque noeud
a exactement 4 fils) conviendra mieux. On définira une classe Matrice
contenant cinq variables :
int n : la profondeur dans l'arbre depuis les feuilles, i.e.
la taille de la matrice est 2n,
Matrice M11, M12, M21, M22 : les 4 sous-matrices correspondants
aux blocs haut-gauche, haut-droit, bas-gauche, bas-droit,
int val : la valeur de l'entrée de la matrice (lorsqu'on est sur une
feuille de l'arbre, i.e.n=0). Pour la partie inversion de matrice, on utilisera
un type float pour val.
On écrira une méthode permettant de passer de la représentation tableau à la représentation
récursive, et inversement.
Tout d'abord, une bonne référence autour de thèmes reliés est le livre
d'Aho, Hopcroft et Ullman [AHU75].
Avant de considérer l'une de ces extensions, il est fortement conseillé de prendre contact.
Des algorithmes plus performant mais plus sophistiqués pour la multiplication de matrices existent.
Ils n'utilisent que O(N2.5) opérations [CW81] ou encore O(N2.38) opérations [CW87].
Hélas leur implémentation est délicate.
Lorsque les matrices sont creuses, il y a des algorithmes meilleurs qui atteignent de meilleurs bornes [YZ04].
Tous ces algorithmes restent valident sur un corps quelconque. On pourra les
adapter pour les corps finis, ou encore sur les rationnels sans avoir à passer par le type
des réels.
Les décompositions standards (LU, LUP...) ont le même coût que la
multiplication rapide. L'algorithme utilise une version récursive par blocs
de la décomposition [AHU75]. Il est alors possible de calculer valeurs propres,
déterminant, et inverse dans la foulée.
La recherche de couplage parfait peut être généralisée à un graphe quelconque
soit en utilisant encore la multiplication de matrices [MVV87],
soit en utilisant la technique de maximisation de flots [MV80].
Enfin, pour la détection de k-cycles,
il y a une amélioration consistant à colorier aléatoirement le graphe avec k couleurs
plutôt qu'en l'orientant. On gagne au niveau de la probabilité de succès puisque k! est alors
remplacé par 2k [AYZ95]. Dans le cas de petits cycles comme le triangle,
ou si le cycle est de taille paire ou encore un multiple de 4, il y a d'autre
stratégies plus rapide [AYZ97].
L'adresse
http://www.enseignement.polytechnique.fr/profs/informatique/Frederic.Magniez/04-05/INF_431/projet.html
contient certains des articles référencés. Si nécessaire, des indications complémentaires y
seront ajoutées au fur et à mesure.
Alfred V. Aho, John E. Hopcroft, and Jeffrey D. Ullman.
The design and analysis of computer algorithms.
Reading, Mass. : Addison-Wesley Pub. Co., 1975 printing, 470 p. CALL
NUMBER: QA76.6 .A36, 1975.
Coppersmith, D. and Winograd, S.
Matrix multiplication via arithmetic progression.
In Proc. 19th Annual ACM Symposium of Theory of
Computing, pages 1–6, 1987.
Rusins Freivalds.
Fast probabilistic algorithms.
In J. Becvár, editor, Mathematical Foundations
of Computer Science 1979, volume 74 of lncs, pages 57–69, Olomouc,
Czechoslovakia, 3–7 September 1979. Springer-Verlag.
S. Micali and V. Vazirani.
An algorithm for maximum matching in general graphs.
In Proceedings of the 21st IEEE Annual Symposium on the
Foundations of Computer Science, pages 17–27, 1980.