\[
\def\CC{\bf C}
\def\QQ{\bf Q}
\def\RR{\bf R}
\def\ZZ{\bf Z}
\def\NN{\bf N}
\]
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
La matrice inverse:
Par pivot de Gauß:
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,
\[\begin{split}\begin{aligned}
M^{-1} =
\begin{pmatrix}
A^{-1}+A^{-1}B\Delta^{-1}CA^{-1} & - A^{-1}B \Delta^{-1}\\
-\Delta^{-1}CA^{-1} & \Delta^{-1}\\
\end{pmatrix}\,.
\end{aligned}\end{split}\]
Voir:
https://en.wikipedia.org/wiki/Block_matrix#Block_matrix_inversion
Exercice
Vérifier que l’on retrouve bien la formule d’inversion des matrices
\(2\times 2\) à partir de la formule d’inversion par blocs.
Vérifier 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\)),
\[d_{2n}
\,\leq\, 2d_n + 14c_n
\,\leq\, 2ban^\omega + 14an^\omega
\,\leq\, \frac{1+\frac 7b}{2^{\omega-1}} \,ba(2n)^\omega\]
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:
\[a_1 = a_0 - \frac{f(a_0)}{f'(a_0)}\]
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
Exercice
Soit \(B(z)\) une série avec un terme constant non nul. Elle admet alors
une unique série inverse \(A(z)\).
Soit \(C(z)\) une série. Vérifier 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?
Proposer 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!
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)\).
Proposition: Si \(c_n=O\left(nf(n)\right)\), alors
\(d_n=O\left(nf(n)\right)\).
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\).
\[d_{2n}
\,\leq\, d_n + 2c_{2n}
\,\leq\, banf(n) + 2a(2n)f(2n)
\,\leq\, (\frac12 + \frac 2b) \,ba(2n)f(2n)\]
On conclue 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 \(\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, proposer 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?
Proposer 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
\[F = GQ + R\]
Récrivons ceci sous la forme:
\[\frac FG = Q + \frac RG\]
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 z^{g-r}\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]\)
:
\[A = A(z) = a_0 + a_1 z + \cdots + a_n z^n\]
\[B = B(z) = b_0 + b_1 z + \cdots + b_n z^m\]
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
Donner 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:
Étape de récurrence
Supposons que \(n=m=2l\), et écrivons
\[A = A_0 + A_1 z^l\]
\[B = B_0 + B_1 z^l\]
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\) :
\[AB = A_0B_0 + ( A_0B_1 + A_1B_0 ) z^l + (A_1B_1)z^{2l}\]
Ou seulement avec trois:
\[AB = A_0B_0 + ( (A_0+A_1)(B_0+B_1) - A_0B_0 - A_1B_1 ) z^l + (A_1B_1)z^{2l}\]
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:
\[f(r)=3f(r-1)=3^rf(0)=3^r\]
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:
\[3^{\lceil \log_2 n \rceil} \leq 3. 3^{\log_2 n} = 3.2^{\log_2 3 . \log_2 n} = 3 n^{\log_2 3}\]
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
\(\NN[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:
\[C(x_0) = A(x_0) B(x_0)\]
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:
\[\begin{split}\begin{aligned}
\Phi: \begin{cases} K[z] &\mapsto (K^n,+,.)\\ P(z) &\mapsto ( P(x_1), \ldots, P(x_n) ) \end{cases}
\end{aligned}\end{split}\]
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:
\[A(z)B(z) = \Phi^{-1} ( \Phi(A) \Phi(B) )\]
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?
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 …
