Ce document dans d'autres formats: feuille de travail, source RST.
Dans ce cours, nous explorons quelques limites de l'algorithme de Gauß en terme de performances, et donnons des approches pour les contourner.
On a vu que l'algorithme de Gauß était de complexité $O(n^3)$, au moins dans son implantation naïve. Vérifions cela expérimentalement avec les petits outils suivants:
load("media/gauss.py")
m = matrice_inversible(3, GF(7)); m # random
gauss(m) # random
Commençons par un corps fini:
import functools
construit_donnee = functools.partial(matrice_inversible, corps=GF(7))
construit_donnee(3) # random
temps(gauss, 10, construit_donnee)
temps(gauss, 20, construit_donnee)
t = [temps(gauss, 2^k, construit_donnee) for k in range(9)]; t
[ t[i+1]/t[i] for i in range(len(t)-1) ]
Un peu plus de valeurs supplémentaires calculées au préalable:
t = [2.8848648071289062e-05, 2.9087066650390625e-05, 4.1961669921875e-05, 9.298324584960938e-05, 0.0003731250762939453, 0.0017020702362060547, 0.010125160217285156, 0.04890704154968262, 0.3750150203704834, 2.7238361835479736, 20.545907974243164, 182.26634407043457, 1334.3144991397858]
[ t[i+1]/t[i] for i in range(len(t)-1) ]
points(enumerate(_))
C'est plausible.
Prenons maintenant le corps des rationnels:
construit_donnee = functools.partial(matrice_inversible, corps=QQ)
t = [temps(gauss, 2^k, construit_donnee) for k in range(8)]; t # random
t = [6.389617919921875e-05, 0.00010395050048828125, 0.000308990478515625, 0.001764059066772461, 0.012479066848754883, 0.09727597236633301, 0.8789999485015869, 9.599533081054688, 127.58634281158447, 2059.1530270576477]
[ t[i+1]/t[i] for i in range(len(t)-1) ]
Bof ...
Avec la matrice de Hilbert:
def hilbert(n):
return matrix(QQ, n, n, lambda i,j: 1/(1+i+j))
hilbert(3)
[temps(gauss, 2^k, hilbert) for k in range(8)] # random
t = [2.193450927734375e-05, 3.4809112548828125e-05, 9.202957153320312e-05, 0.0004980564117431641, 0.003587961196899414, 0.029154062271118164, 0.2275228500366211, 2.2509679794311523, 25.5708429813385, 328.3195171356201]
[ t[i+1]/t[i] for i in range(len(t)-1) ]
Bof!
Prenons un corps de fractions rationnelles:
K = QQ['x'].fraction_field()
construit_donnee = functools.partial(random_matrix, K)
construit_donnee(2)
t = [temps(gauss, n, construit_donnee) for n in range(10)]; t
t = [1.1920928955078125e-05, 1.9073486328125e-05, 0.00018310546875, 0.0007388591766357422, 0.002237081527709961, 0.007543087005615234, 0.021083831787109375, 0.08204507827758789, 0.15540504455566406, 0.9841179847717285, 22.702343940734863]
[ t[i+1]/t[i] for i in range(len(t)-1) ]
Pourquoi est-ce que notre analyse de complexité est si éloignée de la réalité?
Parce que l'on a un mauvais modèle!
On a mesuré la complexité arithmétique de l'algorithme de Gauß, la métrique étant donnée par le nombre d'opérations arithmétiques.
Or, comme l'a constaté toute personne ayant calculé un pivot de Gauß à la main, les coefficients ont tendance à grossir:
%display latex # not tested
gauss(matrice_inversible(10)) # random
def max_coeff(m): return max([c.numer() for row in m.rows() for c in row])
t = [ max_coeff(gauss(matrice_inversible(2^k))) for k in range(7) ]
t
Considérer que le temps nécessaire à une opération arithmétique est constant est donc abusif!
Une meilleure mesure est donc la complexité en bits, puisque les opérations arithmétiques sont de complexité polynomiale en $n$ (en fait en gros $n\log n$) où $n$ est le nombre de bits:
tt = [x.ndigits(2) for x in t]
[float(tt[i+1]/tt[i]) for i in range(len(t)-1)]
Cela suggère expérimentalement que, pour les rationnels, le nombre de bits est borné par un petit polynôme en $n$.
Exemple
Entrée: une matrice $M$ à coefficients entiers
Sortie: son déterminant
C'est un problème typique: on a un résultat qui est relativement petit (un nombre) mais son calcul direct nécessite de manipuler pleins de gros coefficients.
Algorithme modulaire
Algorithme multimodulaire
Intérêt du multimodulaire?
Remarque
$$\operatorname{rang} (M\mod p) \leq \operatorname{rang}(M)$$La clef des algorithmes précédents est que l'on avait un morphisme du corps de base dans un/plusieurs corps où l'arithmétique était plus efficace, avec la possibilité d'inverser localement ce morphisme à la fin.
La même technique s'adapte à chaque fois que l'on a une explosion des coefficients intermédiaires, et un morphisme dans un (ou plusieurs) corps où l'on maîtrise la croissance des coefficients intermédiaires.
Exemple de problème
Entrée: une matrice $M$ à coefficients polynomiaux
Sortie: son déterminant
Algorithme
Exercice
Donner une borne de complexité pour le calcul du polynôme caractéristique d'une matrice dans $GF(p)$.
Exemple de problème
Entrée: une matrice $M$ carrée inversible à coefficients entiers,
Sortie: l'inverse de $M$
On considère le morphisme partiel $\phi$ de $\QQ$ dans $GF(p^k)$ :
Algorithme
Raffiner itérativement cette solution:
$k$ double à chaque itération!!!
Mise en contexte: on a écrit notre matrice $M$ comme une série:
$$M = M_0 + M_1 p + M_2 p^2 + \cdots$$où chaque $M_i$ est essentiellement une matrice mod $p$, et on a utilisé la technique classique de l'itération de Newton pour calculer une solution de plus en plus précise de l'équation $MN=1$.
Considérons une matrice creuse:
M = random_matrix(GF(7), 19, sparse=True, density=1/3)
M
Et appliquons un pivot de Gauß partiel:
gauss(M,10)
Regardons à plus grande échelle:
M = random_matrix(GF(7), 200, sparse=True, density=1/10)
len(M.nonzero_positions())
len(gauss(M, 50).nonzero_positions())
Problème
Beaucoup de matrices apparaissant dans les problèmes pratiques sont structurées:
L'algorithme de Gauß ne préserve pas ces structures!
Et pourtant:
M = random_matrix(GF(7), 10000, sparse=True, density=3/10000)
M.rank()
Comment cela marche-t-il???
Problème
Calculer le polynôme minimal d'une matrice
Remarque
Soit $P$ le polynôme minimal d'une matrice carrée $M$.
Soient $U$ et $V$ deux vecteurs.
Alors la suite de nombre $u_k = U M^k V$ satisfait une relation de récurence donnée par les coefficients de $P$.
Rappel: Algorithme de Berlekamp-Massey
L'algorithme de Berlekamp-Massey permet, étant donné une suite $s_{1},\dots,s_{n}$ d'éléments d'un corps de trouver la plus petite relation de récurrence satisfaite par cette suite. Les coefficients de cette relation de récurrence sont traditionnellement encodés sous la forme d'un polynôme. Encore une conséquence d'Euclide étendu ...
Voir TP pour les détails.
Algorithme de Wiedemann
Remarques
Exercice
Supposer que le polynôme minimal de $M$ soit $X^3-2X+1$.
Déterminer l'inverse de $M$.
Voir TP.
http://en.wikipedia.org/wiki/Newton%27s_identities
Exercice: Berlekamp-Massey
L'algorithme de Berlekamp-Massey permet, étant donné une suite $s_{1},\dots,s_{n}$ d'éléments d'un corps de trouver la plus petite relation de récurrence satisfaite par cette suite. Les coefficients de cette relation de récurrence sont traditionnellement encodés sous la forme d'un polynôme.
Cette algorithme est implanté dans Sage par la fonction berlekamp_masse. Vous pouvez au choix faire quelques essais avec cette fonction et passer à l'exercice suivant ou ...
Implanter l'algorithme de Berlekamp-Massey, soit en vous servant de [Massey.1969]: Shift-register synthesis and BCH Decoding James L. Massey, IEEE transactions on information theory, 1969, ou via l'algorithme d'Euclide étendu décrit dans le texte sur Wiedemann ou, avec plus de détails dans le texte sur le code de Goppa.
Exercice: Polynome minimal et Wiedemann sur un exemple
diagonalizable
; on pourra tirer les multiplicités au hasard avec IntegerVectors.minimal_polynomial
.Exercice: Implantation de l'algorithme de Wiedemann
endomorphisme
qui prend une matrice $M$ en argument, et renvoie l'endomorphisme correspondant, c'est-à-dire la fonction qui à un vecteur $v$ associe le vecteur $M\times v$.wiedemann(f, V)
qui prend un endomorphisme $f$ et l'espace sur lequel il agit, et calcule son polynôme minimal en utilisant l'algorithme de Wiedemann.Exercice: Calcul du rang par préconditionnement par produit de matrices diagonales.
On s'intéresse aux graphes simples (pas de boucles, pas d'arêtes multiples, pas d'orientation, ...) à isomorphie près (les sommets ne portent pas d'étiquette permettant de les distinguer). Une forêt est un graphe acyclique; un arbre est une forêt connexe.
Exercice
Fabriquer la liste de toutes les forêts à $6$ sommets et $4$ arêtes, à isomorphie près. Que constatez-vous?
Indication: essayer:
for g in graphs(5): show(g)
puis utiliser les options property
et size
et la méthode is_forest
.
Fixons un entier $n$. On va considérer la matrice $T_n$ dont
Exercice
Écrire une fonction qui construit la matrice $T_n$.
Indications:
G = Graph([[1,2]])
{G: 1}
G = G.copy(immutable=True)
{G: 1}
copy
et delete_edge
pour obtenir $f$ par suppression d'une arête de $t$. Puis utiliser la méthode canonical_label
pour mettre $f$ sous forme canonique à isomorphie près.Exercice
Explorer les propriétés de la matrice $T_n$.
Pour les curieux, voir arXiv:0912.2619, p. 21.