Produits rapides#
Nicolas M. Thiéry <Nicolas.Thiery at universite-paris-saclay.fr>
Motivation: «Tout se ramène aux produits»#
Inversion de matrices#
Exercice: matrices \(2\times 2\) génériques
Soit \(M=\begin{pmatrix}a&b\\c&d\end{pmatrix}\).
Calculer \(M^{-1}\) en utilisant les cofacteurs.
Calculer \(M^{-1}\) par pivot de Gauß.
Refaire le même calcul avec une matrice par blocs \(M=\begin{pmatrix}A&B\\C&D\end{pmatrix}\)!
Solution
%display latex
a,b,c,d = QQ['a,b,c,d'].fraction_field().gens()
M = matrix([[a,b],[c,d]]); M
La matrice inverse:
M^-1
Par pivot de Gauß:
I2 = matrix(2,2,1); I2
M = M.augment(I2, subdivide=True); M
M[1] = M[1] - c/a *M[0]; M
M[1] = M[1]/M[1,1]; M
M[0] = M[0] - b * M[1]; M
M[0] = M[0]/a; M
Théorème: formule d’inversion de matrice par blocs
Soit \(M=\begin{pmatrix}A&B\\C&D\end{pmatrix}\) une matrice par blocs, où \(A\) et \(D\) sont carrées. On suppose de plus que \(A\) et \(\Delta=(D-CA^{-1}B)\) sont inversibles. Alors,
Exercice
Vérifiez que l’on retrouve bien la formule d’inversion des matrices \(2\times 2\) à partir de la formule d’inversion par blocs.
Vérifiez que l’on obtient bien l’inverse de \(M\).
Théorème: «pour les matrices, l’inversion ne coûte pas plus que la multiplication»
Soient respectivement \(c_n\) et \(d_n\) les complexités de la multiplication et de l’inversion de matrices de taille \(n\).
Si \(c_n=O(n^\omega)\) avec \(\omega\geq 2\), alors \(d_n=O(n^\omega)\).
Démonstration (simplifiée):
Soit \(a\) tel que \(c_n\leq a n^\omega\), pour tout \(n\). On va chercher \(b\) tel que \(d_n \leq bc_n=ba n^\omega\), pour tout \(n\).
Supposons que \(n\) et \(b\) sont tels que \(d_n \leq ba n^\omega\).
Alors, en utilisant la formule par blocs ci-dessus, et en comptant les additions avec les multiplications (\(\omega \geq 2\)),
Donc, à condition de prendre \(b\) suffisamment grand, \(d_{2n} \leq ba(2n)^\omega\), comme voulu.
Cela donne par récurrence la complexité voulue pour \(n\) une puissance de \(2\). Quitte à gérer proprement les parties entières, le même argument marche pour tout \(n\).
Méthode de Newton#
Approximation numérique de solutions de \(f(a) = 0\)#
Exercice
Soit \(f(x)\) une fonction suffisamment gentille dont on recherche une racine \(a\).
On suppose que l’on dispose d’une approximation \(a_0\) de \(a\), et on pose:
Calculer \(f(a)\) par développement de Taylor de \(f\) en \(a_0\).
Qu’en déduire sur \(a_1-a\) par rapport à \(a_0-a\)?
Quelle conclusion peut-on en tirer? Sous quelles hypothèses?
Pour les détails, voir l’article de la Wikipedia.
Inversion de séries formelles#
Exercice
Soit \(B(z)\) une série formelle avec un terme constant non nul. Elle admet alors une unique série inverse \(A(z)\).
Soit \(C(z)\) une série. Vérifiez que:
\[C(z) = A(z) + O(z^k) \Longleftrightarrow C(z)B(z) = 1 + O(z^k)\]On suppose que l’on dispose d’une approximation \(A_0(z)\) de l’inverse \(A(z)\) de \(B(z)\) :
\[A_0(z)B(z)=1+O(z^k)\]et on pose:
\[A_1(z)=A_0(z)(2-A_0(z)B(z))\]Que peut on dire de cette nouvelle approximation?
Proposez un analogue pour l’inversion des séries du théorème sur l’inversion des matrices vu précédemment.
On verra en TP que l’expression ci-dessus pour \(A_1(z)\) peut être obtenue par itération de Newton.
Solution
Pour \(\Longrightarrow\) :
\[C(z)B(z) - 1 = (C(z)-A(z)) B(z) = O(z^k) B(z) = O(z^k)\]Pour \(\Longleftarrow\) :
\[C(z)- A(z) = (C(z)B(z)-1) A(z) = O(z^k) A(z) = O(z^k)\]Posons maintenant \(C(z)\) tel que \(A_0(z)B(z) = 1 + C(z)\), et calculons \(A_1(z)B(z)\) :
\[A_1(z)B(z) = A_0(z)(2-A_0(z)B(z))B(z) = (1+C(z)) (1-C(z)) = 1-C(z)^2 = 1 + O(z^{2k})\]La précision double!
Proposition
Notons \(c_n\) la complexité de la multiplication et \(d_n\) la complexité de la division de séries. Soit \(f(n)\geq 0\) croissante; par exemple \(f(n)=n^\omega\) avec \(\omega>0\) ou \(f(n)=\log(n)\).
Si \(c_n=O\left(nf(n)\right)\), alors \(d_n=O\left(nf(n)\right)\).
Démonstration
Soit \(a\) tel que \(c_n\leq a nf(n)\), pour tout \(n\). On va chercher \(b\) tel que \(d_n \leq ba nf(n)\), pour tout \(n\).
On conclut comme précédemment.
Approximation en série d’une équation implicite \(F(A(z)) = 0\)#
Problème: calcul de la racine d’une série
Soit \(B(z)\) une série. On souhaite calculer sa racine, c’est-à-dire la série \(A(z)\) telle que \(A(z)^2-B(z)=0\).
Comment procéder?
Exercice
Soit \(F(X)\) un polynôme à coefficients dans \(\mathbb{QQ}[[z]]\). Par exemple: \(F(X)=X^2 - B(z)\).
On cherche une série \(A(z)\) telle que \(F(A(z))=0\).
On suppose que l’on dispose d’une approximation \(A_0(z)\) de \(A(z)\).
En vous inspirant de la méthode de Newton usuelle, proposez une meilleure approximation \(A_1(z)\) de \(A(z)\).
Quelle est la vitesse de convergence?
Quelles opérations sont nécessaires lors d’une itération?
Proposez un analogue pour la résolution des équations implicites du théorème sur l’inversion de matrices vu précédemment.
Exercice
En déduire un algorithme pour calculer la racine carrée d’une série \(B(z)\). Que faut-il comme hypothèse sur \(B(z)\)?
Que se passe-t’il si l’on essaye de calculer l’inverse d’une série de cette manière?
Remarque
On peut en fait traiter de manière similaire des équations différentielles linéaires \(F(A(z))=0\).
Application: connaissant \(B(z)\), calculer \(A(z)=\exp(B(z))\)
Indication: Prendre \(F(A(z)) = A(z)'-B(z)'A(z)\)
Division Euclidienne de polynômes#
Soit \(F=F(z)\) et \(G(z)\) deux polynômes. On veut déterminer \(Q\) et \(R\) avec \(\deg R < \deg G\) tels que
Récrivons ceci sous la forme:
Que se passe-t’il si \(z\) est «grand»?
Idée:
Envoyer l’infini sur zéro en prenant la réciproque des polynômes concernés: \(\overline P = z^{\deg P} P(\frac1z)\)
Notant \(f,g,q,r\) les degrés des polynômes, et remarquant que \(f=g+q\), \(r<g\) et \(\overline G\) a une constante non nulle, on obtient: $\( \frac{\overline F}{\overline G} = \overline Q + z^{q+g-r}\frac{\overline R}{\overline G}=Q + O(z^{q+1})\frac{\overline R}{\overline G}\)$.
Obtenir \(\overline Q\) en calculant les premiers termes de \(\frac{\overline F}{\overline G}\). C’est une inversion de série: itération de Newton, …
Calculer \(R=F-GQ\)
Pour les détails, voir page 83 de [AECF]
Produits rapides#
Produits rapides de polynômes#
Dans la suite, on considère un anneau \(K\) et deux polynômes dans \(K[z]\) :
L’objectif est de calculer les coefficients \(c_k\) du polynôme \(C(z) = A(z)B(z)\).
Algorithme naïf#
On se contente d’utiliser la formule \(c_k = \sum_{i+j=k} a_i b_j\).
Exercice
Quelle est la complexité du calcul du produit des polynômes \(A(z)\) et \(B(z)\) par l’algorithme naïf?
Karatsuba#
Exercice
Donnez des formules pour calculer les coefficients du polynôme \((a_0+a_1z)(b_0+b_1z)\) en fonction de \(a_0,a_1,b_0,b_1\) utilisant un nombre minimal de produits.
Nous allons maintenant appliquer les deux principes suivants:
«Si vous avez une bonne idée, appliquez la par récurrence, vous obtiendrez une meilleure idée.»
Diviser pour régner!
Étape de récurrence
Supposons que \(n=m=2l\), et écrivons
où \(A_0=A_0(z)\), \(A_1=A_1(z)\), \(B_0=B_0(z)\) et \(B_1=B_1(z)\) sont de degré \(\leq l\).
On peut calculer \(AB\) en calculant récursivement quatre produits de polynômes de degré \(l\) :
Ou seulement avec trois:
L’algorithme de Karatsuba consiste à calculer le produits de polynômes de degré \(2^r\) en appliquant récursivement l’étape précédente.
Complexité
L’algorithme de multiplication de Karatsuba est de complexité \(O(n^{\log_2(3)})\approx O(n^{1.59})\).
Démonstration
On suppose d’abord que \(n=2^r\), et on ne compte que le nombre \(f(r)\) de multiplications requises dans \(K\). Clairement:
Pour calculer le produit de deux polynômes de degré \(n\), on les complète en polynômes de degré \(2^{\lceil \log_2(n)\rceil}\). Le nombre de multiplications dans \(K\) est alors borné par:
Il est clair que le nombre d’additions est négligeable (de l’ordre de \(O(4n\log_2 n)\)).
En pratique: implantation
L’algorithme de Karatsuba étant plus compliqué, en particulier à cause de la récursion, est moins performant en petit degré que l’algorithme naïf. Aussi les implantations utilisent l’étape de récurrence en haut degré, et basculent sur un produit naïf en deçà d’un certain seuil.
Ce seuil est déterminé expérimentalement par bancs d’essais. Dans certains cas la détermination du seuil optimal pour une architecture donnée est effectuée automatiquement à la compilation.
C’est un principe très général. On l’avait déjà vu avec les tris, et on le retrouve par exemple en algèbre linéaire avec la bibliothèque ATLAS (Automatically Tuned Linear Algebra Software)
En pratique: usage
L’algorithme de Karatsuba requiert des soustractions:
Il ne s’applique pas aux polynômes sur des semi-anneaux (par exemple \(\mathbb{N}[x]\), algèbre tropicale, …)
Il peut poser des problèmes de stabilité numérique en calcul approché (flottants, …)
Produit par évaluation#
Remarque stupide
Si \(x_0\) est un élément de \(K\), et \(C(z) = A(z)B(z)\) alors:
Corollaire
Soient \(x_1,\dots,x_n\) des éléments de \(K\) et munissons \(K^n\) de l’addition et de la multiplication point à point.
L’application d’évaluation:
est un morphisme d’algèbre.
C’est même un isomorphisme si on se restreint à l’ensemble \(K[z]_n\) des polynômes de degré \(<n\).
Le produit dans \((K^n,+,.)\) est de complexité \(n\). Donc il est tentant d’utiliser cet isomorphisme pour calculer les produits:
Problème
Rentable si le calcul de \(\Phi\) (évaluation) et de \(\Phi^{-1}\) (par ex. interpolation) est peu coûteux. Pour des points quelconques, c’est au moins du \(O(n^2)\).
Comment choisir de bons points d’évaluation?
Produit par transformée de Fourier Discrète#
Transformée de Fourier Discrète#
Proposition
Supposons que l’anneau \(K\) contienne une racine primitive \(\omega\) de l’unité. Alors le morphisme d’algèbre:
induit un isomorphisme d’algèbre de \(K[z] / (z^n-1)\).
Démonstration
Regarder le noyau + dimension.
Remarque
On retrouve la même algèbre que dans les codes cycliques; entre autres, la multiplication par \(z\) donne une action du groupe cyclique \(C_n\).
Exercice
\(DFT_\omega\) est une application linéaire. Donnez sa matrice.
Donnez la matrice inverse.
Indication: \(\sum_{k=0}^{n-1} \omega^{ik} = \begin{cases}n&\text{ si } i\equiv 0[n]\\0&\text{ sinon}\end{cases}\)
Proposition
La transformée de Fourier discrète inverse est encore une transformée de Fourier discrète:
Remarque: lien avec la théorie des représentations
La matrice de \(DFT_\omega\) est aussi la table des caractères du groupe cyclique \(C_n\). Le fait qu’elle soit unitaire à un scalaire près est un cas particulier d’une proposition générale sur les tables de caractères. L’espace \(K[z]/(z^n-1)\) se décompose en \(n\) modules simples de dimension \(1\) et la transformation \(DFT_\omega\) correspond à la décomposition d’un polynôme dans ces modules simples.
Il existe des notions de transformées de Fourier discrètes pour d’autres groupes.
Il reste à calculer efficacement la transformée de Fourier discrète.
Transformée de Fourier rapide (FFT: Fast Fourier Transform)#
Diviser pour régner
Supposons que \(P\) soit un polynôme de degré au plus \(n=2k\).
Noter que \(z^{2k} - 1 = (z^k-1) (z^k+1)\).
Du coup, la moitié des racines \(2k\)-ièmes sont des racines \(k\)-èmes de l’unité, racines de \(z^k-1=0\). On peut donc utiliser la transformée de Fourier discrète pour évaluer \(P(z)\) dessus. Plus précisément, on calcule
(ce calcul est léger!) et on utilise \(DFT_{\omega^2}(P_+(z))\) pour retrouver l’évaluation de \(P(z)\) aux racines \(k\)-èmes de l’unité.
L’autre moitié des racines \(2k\)-ièmes sont les racines \(k\)-ème de l’unité décalées par un facteur \(\omega\), racines de \(z^k+1\). On calcule alors
et on peut donc utiliser \(DFT_{\omega^2}(P_-(\omega z))\).
Algorithme de multiplication par FFT
On considère une racine \(2^k\)-ème de l’unité, et on applique récursivement l’idée précédente.
Complexité: \(O(n\log n)\), comme pour les tris.
Problème
Et s’il n’y a pas de racine primitive de l’unité dans \(K\)?
On la rajoute!
Exemple: les corps cyclotomiques obtenus par extension algébrique de \(\mathbb{Q}\) par un polynôme cyclotomique:
K = CyclotomicField(6)
omega = K.gen()
omega^6
Souci: ces corps cyclotomiques nécessitent de calculer dans des extensions de corps de haut degré; donc un bon produit; cela pourrait boucler!
Algorithme de Schönhage et Strassen: \(O(n\log n\log\log n)\)
Autre souci: on a divisé par \(n=2^k\); ce n’est pas forcément possible, par exemple en caractéristique \(2\)!
Produits rapides d’entiers
Même principe que pour les polynômes; juste plus technique à cause de la gestion des retenues. On retrouve le produit par Karatsuba, par FFT, …
Ce que l’on a remarqué pour les séries s’applique aux calculs sur les nombres réels à précision arbitraire.
Produits rapides de matrices
Algorithme de Strassen
Même principe que Karatsuba!
Pour multiplier deux matrices \(2\times 2\), il existe des formules n’utilisant que 7 produits au lieu de 8.
On découpe les matrices de taille \(2^k\) en \(4\) blocs de taille \(2^{k-1}\) et on utilise les formules ci-dessus récursivement.
Complexité: \(O(n^{\log_2 7}) \approx O(n^{2,8})\)
Améliorations: pour un \(k\) donné (ci-dessus \(k=2\)), chercher systématiquement l’algorithme optimal. Puis l’appliquer récursivement en découpant les matrices en \(k\times k\) blocs.
Algorithmes de Coppersmith-Winograd et suivants
Complexité: \(O(n^{2.3755\cdots})\) (1990), \(O(n^{2.374})\) (2010), \(O(n^{2,3728\cdots})\) (2014), …
Inutilisables en pratique …
Résumé#
«Tout se ramène aux produits»
calcul d’inverse
division Euclidienne
résolution d’équation (différentielle)
Algorithmes récursifs (diviser pour régner):
Structures récursives: matrices et polynômes par blocs
Itération de Newton
Évaluation, FFT