\[ \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 u-psud.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}\).

  1. Calculer \(M^{-1}\) en utilisant les cofacteurs.

  2. Calculer \(M^{-1}\) par pivot de Gauß.

  3. Refaire le même calcul avec une matrice par blocs \(M=\begin{pmatrix}A&B\\C&D\end{pmatrix}\)!

Solution

La matrice inverse:

%display latex
a,b,c,d = QQ['a,b,c,d'].fraction_field().gens()
M = matrix([[a,b],[c,d]]); M
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,

\[\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

  1. Vérifier que l’on retrouve bien la formule d’inversion des matrices \(2\times 2\) à partir de la formule d’inversion par blocs.

  2. 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)}\]
  1. Calculer \(f(a)\) par développement de Taylor de \(f\) en \(a_0\).

  2. Qu’en déduire sur \(a_1-a\) par rapport à \(a_0-a\)?

  3. 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)\).

  1. 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)\]
  2. 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?

  3. Proposer un analogue pour l’inversion des séries du théorème sur l’inversion de matrices ci-dessus.

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)\).

  1. En vous inspirant de la méthode de Newton usuelle, proposer une meilleure approximation \(A_1(z)\) de \(A(z)\).

  2. Quelle est la vitesse de convergence?

  3. Quelles opérations sont nécessaires lors d’une itération?

  4. Proposer un analogue pour la résolution des équations implicites du théorème sur l’inversion de matrices ci-dessus.

Exercice

  1. 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)\)?

  2. 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\), \(q<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:

  • «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

\[A = A_0 + A_1 z^l\]
\[B = B_0 + B_1 z^l\]

\(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?

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:

\[\begin{split}\begin{aligned} DFT_\omega: \begin{cases} K[z] &\mapsto (K^n,+,.)\\ P(z) &\mapsto ( P(1), P(w), \ldots, P(w^{n-1}) ) \end{cases} \end{aligned}\end{split}\]

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

  1. \(DFT_\omega\) est une application linéaire. Donner sa matrice.

  2. Donner 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:

\[DFT_\omega^{-1} = \frac 1n DFT_{\omega^{-1}}\]

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

\[P_+(z) = P(z) [ z^k - 1 ]\]

(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

\[P_-(z) = P(z) [ z^k + 1 ]\]

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 \(\QQ\) 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

    • FFT