**Théorème 8.2** — La multiplication matricielle naïve n'est pas optimale. Soit $\mathbb{K}$ un corps. On peut multiplier deux matrices de $M_n(\mathbb{K})$ en

a. $O(n^3)$ opérations dans $\mathbb{K}$, par l'algorithme naïf;

b. $n^2[\frac{n}{2}] + 2n[\frac{n}{2}] \approx \frac{1}{2}n^3 + n^2$ multiplications dans $\mathbb{K}$;

c. $n^2[\frac{n}{2}] + (2n-1)[\frac{n}{2}] \approx \frac{1}{2}n^3 + n^2 - \frac{n}{2}$ multiplications dans $\mathbb{K}$ si la caractéristique de $\mathbb{K}$ est différente de 2 et si la division par 2 est gratuite;

d. $O(n^{\log_2 7}) \approx O(n^{2,81})$ opérations dans $\mathbb{K}$.

Ces résultats restent valables sur un anneau commutatif $\mathbb{A}$; pour (a) et (d) l'hypothèse de commutativité peut même être supprimée. La partie (b) est due à Winograd et fut historiquement la première amélioration de l'algorithme naïf. Elle réduit de moitié le nombre de multiplications, en augmentant en revanche le nombre d'additions. Cela constitue déjà un progrès pour une large classe d'anneaux $\mathbb{A}$ dans lesquels la multiplication est plus coûteuse que l'addition $^2$. Par exemple, la technique du « scindage binaire » qui sera utilisée au Chapitre 15 pour le calcul d'un terme d'une récurrence, ou encore pour le calcul du développement décimal de $e$, requiert de multiplier des matrices contenant de gros nombres entiers.

L'amélioration (c) a été obtenue par Waksman, et implique qu'on peut effectuer le produit de deux matrices $2 \times 2$ en 7 multiplications. L'inconvénient commun des algorithmes de Winograd et de Waksman est l'utilisation de la commutativité de $\mathbb{K}$, qui devient un obstacle à une application récursive, ne permettant pas d'améliorer (la borne supérieure 3 sur) l'exposant $\omega_{\text{mul}}$. Beaucoup pensaient que tout résultat de type (b) et (c) serait optimal, au sens que $\frac{1}{2}n^3$ multiplications dans $\mathbb{K}$ seraient nécessaires pour le produit dans $M_n(\mathbb{K})$. Ceci fut contredit par la découverte par Strassen de (d); ce résultat découle d'un fait d'une simplicité déconcertante, mais d'une portée considérable : deux matrices $2 \times 2$ peuvent être multipliées en 7 multiplications même si $\mathbb{K}$ n'est pas commutatif.

La non-optimalité de la multiplication matricielle naïve justifie l'introduction de la définition suivante.

**Définition 8.2** — Fonction de multiplication matricielle. On dit que l'application $\text{MM}: \mathbb{N} \to \mathbb{N}$ est une fonction de multiplication matricielle (pour un corps $\mathbb{K}$) si :

- on peut multiplier deux matrices arbitraires de $M_n(\mathbb{K})$ en au plus $\text{MM}(n)$ opérations dans $\mathbb{K}$;

- $\text{MM}$ est croissante et vérifie l'inégalité $\text{MM}(n/2) \le \text{MM}(n)/4$ pour $n \in \mathbb{N}$.

Ainsi, le Théorème 8.2 montre que les fonctions $n \mapsto n^3$ ou $n \mapsto n^{\log_2 7}$ sont (à des constantes près) des fonctions de multiplication matricielles. Comme dans le cas de la fonction de multiplication polynomiale $M$ introduite au Chapitre 2, l'intérêt de cette notion est qu'elle permet d'énoncer des résultats de complexité indépendants du choix de l'algorithme utilisé pour multiplier les matrices.

La seconde condition de la Définition 8.2 établit la cohérence avec la Définition 8.1 : si $\text{MM}(n) = c \cdot n^\theta$, pour une certaine constante $c > 0$, alors $\theta$ doit être supérieur ou égal

<footnote n="2">2. C'est par exemple le cas de $\mathbb{Z}$ et de $\mathbb{Z}[X]$, mais pas de $\mathbb{Q}$ ou de $\mathbb{Q}(X)$.</footnote>
à 2. Cette condition permettra de montrer que dans des algorithmes de type « diviser pour régner » pour les matrices, le coût total est essentiellement celui du dernier appel récursif.

Le reste de cette section est dédié à la preuve du Théorème 8.2.

## Multiplication naïve

L'algorithme « cubique » de multiplication prend en entrée deux matrices $A = (a_{i,j})_{i,j=1}^n$ et $X = (x_{i,j})_{i,j=1}^n$ de $M_n(\mathbb{K})$, et calcule leur produit $R = AX$ en utilisant la formule

$$r_{i,k} = \sum_{j=1}^{n} a_{i,j} x_{j,k}.$$

Le coût nécessaire pour obtenir un élément de $R$ est donc de $n$ multiplications et $n-1$ additions, de sorte que la complexité totale est en $O(n^3)$.

## Algorithme de Winograd

Il est possible de diviser essentiellement par deux le nombre de multiplications scalaires suffisantes pour le calcul du produit AX. Ceci montre au passage que la multiplication matricielle naïve n'est pas optimale. Supposons pour simplifier que $n$ est pair, $n = 2k$. L'algorithme de Winograd repose sur l'identité suivante, héritière de celle de Karatsuba : si $\ell$ est un vecteur ligne $(a_1 \cdots a_n)$ et $c$ est un vecteur colonne $^t(x_1 \cdots x_n)$, alors le produit scalaire $(\ell | c) = \sum_i a_i x_i$ s'écrit

$$ (\ell|c) = (a_1+x_2)(a_2+x_1) + \cdots + (a_{2k-1}+x_{2k})(a_{2k}+x_{2k-1}) - \sigma(\ell) - \sigma(c), \quad (8.1) $$

où $\sigma(\ell) = a_1a_2 + \cdots + a_{2k-1}a_2k$ et $\sigma(c) = x_1x_2 + \cdots + x_{2k-1}x_2k$. Notons que $\sigma(\ell)$ se calcule en $k = n/2$ multiplications et $k-1 = n/2-1$ additions.

Si maintenant $\ell_1, \dots, \ell_n$ sont les lignes de $A$ et $c_1, \dots, c_n$ sont les colonnes de $X$, l'élément $(i,j)$ du produit AX vaut $(\ell_i|c_j)$. L'idée est alors de précalculer les quantités $\sigma(\ell_i)$ et $\sigma(c_i)$ pour $1 \le i \le n$; ce précalcul demande $n(\frac{n}{2} + \frac{n}{2}) = n^2$ multiplications et $n(\frac{n}{2} - 1 + \frac{n}{2} - 1) = n^2 - 2n$ additions. Une fois déterminés les $\sigma(\ell_i)$ et $\sigma(c_i)$, on peut calculer l'ensemble des éléments $r_{i,j}$ de $R = AX$ grâce à la formule (8.1) en $n^2 \cdot \frac{n}{2} = \frac{1}{2}n^3$ multiplications et $n^2 \cdot (n + (\frac{n}{2} - 1) + 2) = \frac{3}{2}n^3 + n^2$ additions. Au total, le coût de la multiplication de l'algorithme de Winograd est de $\frac{1}{2}n^3 + n^2$ multiplications et $\frac{3}{2}n^3 + 2n^2 - 2n$ additions. Ainsi, on a essentiellement transformé $n^3/2$ multiplications en additions, ce qui est appréciable.

## Algorithme de Waksman

Pour effectuer le produit de matrices $2 \times 2$

$$ R = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \times \begin{pmatrix} x & y \\ z & t \end{pmatrix}, $$

l'algorithme de Winograd ci-dessus utilise 8 multiplications, en écrivant

$$ R = \begin{pmatrix} (a+z)(b+x) - ab - zx & (a+t)(b+y) - ab - ty \\ (c+z)(d+x) - cd - zx & (c+t)(d+y) - cd - ty \end{pmatrix}. $$
Il n'apporte donc aucun gain par rapport à l'algorithme naïf dans cette taille.

L'algorithme suivant, dû à Waksman, peut être vu comme une amélioration de l'algorithme de Winograd. Il effectue le produit de matrices $2 \times 2$ en 7 opérations, dès lors que dans le corps de base $K$ l'élément 2 est inversible et que la division par 2 y est gratuite. L'idée est d'écrire la matrice $R$ comme

$$R = \frac{1}{2} \begin{pmatrix} \alpha_1 - \alpha_2 & \beta_1 - \beta_2 \\ \gamma_1 - \gamma_2 & \delta_1 - \delta_2 \end{pmatrix},$$

où

$$\begin{align*} \alpha_1 &= (a+z)(b+x), & \alpha_2 &= (a-z)(b-x), & \beta_1 &= (a+t)(b+y), & \beta_2 &= (a-t)(b-y), \\ \gamma_1 &= (c+z)(d+x), & \gamma_2 &= (c-z)(d-x), & \delta_1 &= (c+t)(d+y), & \delta_2 &= (c-t)(d-y), \end{align*}$$

et de faire observer que, grâce à l'identité

$$(\gamma_1 + \gamma_2) + (\beta_1 + \beta_2) = (\alpha_1 + \alpha_2) + (\delta_1 + \delta_2) = 2(ab + cd + xz + ty),$$

l'élément $\delta_2$ est une combinaison linéaire des sept autres produits $\alpha_1, \alpha_2, \dots, \delta_1$, et donc n'a pas besoin d'être obtenu lui aussi par un produit.

L'algorithme qui s'en déduit, tout comme celui de Winograd, ne peut être utilisé de manière récursive : la preuve des formules sous-jacentes repose sur le fait que les éléments du corps de base $K$ commutent $(bz = zb, bt = tb, \dots)$. Lors des appels récursifs, les objets à multiplier seraient eux-mêmes des matrices, pour lesquelles la loi de commutation n'est plus vérifiée. On dit donc que les algorithmes de Winograd et de Waksman sont des algorithmes *commutatifs*.

**Exercice 8.1** Finir la preuve des parties (b) et (c) du Théorème 8.2, en généralisant les arguments précédents à des matrices de taille $n$, où $n \in \mathbb{N}$ n'est pas nécessairement pair. ■

### Algorithme de Strassen

L'algorithme de Strassen fut le premier à descendre en-dessous de $O(n^3)$. Tout comme l'algorithme de Karatsuba pour les polynômes, il repose sur le gain d'une multiplication pour les matrices $2 \times 2$, qui devient un gain dans l'exposant lors d'une application récursive. La non-commutativité du produit matriciel se traduit en revenue par une contrainte supplémentaire : afin de se prêter à un emploi récursif, le nouvel algorithme doit, à la différence de l'algorithme de Waksman, satisfaire à une contrainte de non-commutativité.

Soient donc $A$ et $X$ des matrices de taille 2 à multiplier par un algorithme non commutatif. Intuitivement, cela se traduit par le fait que lors des multiplications utilisées par un tel algorithme, les éléments $a,b,c,d$ de $A$ doivent rester à gauche, et les éléments $x,y,z,t$ de $X$ à droite (cela n'est pas le cas pour les algorithmes de Winograd et de Waksman, qui *mélangent* les éléments de $A$ et $X$).

L'algorithme naïf est bien non commutatif, mais demande d'effectuer 8 multiplications. L'algorithme de Strassen que nous allons présenter revient à

— calculer des combinaisons linéaires des éléments de $A$ et des éléments de $X$,
**Entrée** Deux matrices $A, X \in \mathcal{M}_n(\mathbb{K})$, avec $n = 2^k$.

**Sortie** Le produit AX.

1. Si $n = 1$, renvoyer AX.

2. Décomposer $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$ et $X = \begin{pmatrix} x & y \\ z & t \end{pmatrix}$, avec $a, b, c, d$ et $x, y, z, t$ dans $\mathcal{M}_{n/2}(\mathbb{K})$.

3. Calculer récursivement les produits

$$q_1 = a(x + z),$$

$$q_2 = d(y + t),$$

$$q_3 = (d - a)(z - y),$$

$$q_4 = (b - d)(z + t)$$

$$q_5 = (b - a)z,$$

$$q_6 = (c - a)(x + y), \quad q_7 = (c - d)y.$$

4. Calculer les sommes

$$r_{1,1} = q_1 + q_5,$$

$$r_{1,2} = q_2 + q_3 + q_4 - q_5,$$

$$r_{2,1} = q_1 + q_3 + q_6 - q_7,$$

$$r_{2,2} = q_2 + q_7.$$

5. Renvoyer $\begin{pmatrix} r_{1,1} & r_{1,2} \\ r_{2,1} & r_{2,2} \end{pmatrix}$.

**Algorithme 8.1 – Algorithme de Strassen pour multiplier deux matrices.**

— effectuer 7 produits de ces combinaisons linéaires,

— recombiner les résultats pour obtenir les quatre éléments du produit AX.

Explicitement, les formules à appliquer sont données dans les étapes (3) et (4) de l'Algorithme 8.1. Effectuons maintenant le pas récursif. Quitte à rajouter des colonnes et des lignes nulles, on peut supposer que les matrices A et X sont de taille $n = 2^k$, avec $k \in \mathbb{N}$. On procède à une découpe en blocs de A et de X :

$$A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}, \quad X = \begin{pmatrix} x & y \\ z & t \end{pmatrix}.$$

Dans cette écriture $a, b, c, d$ et $x, y, z, t$ sont donc des matrices carrées de taille $n/2$.

Le point-clé est le suivant : grâce à la non-commutativité, les formules données dans l'Algorithme 8.1 pour le cas $2 \times 2$ s'appliquent toujours si $a, b, c, d$ et $x, y, z, t$ sont des matrices. Ainsi, elles permettent de ramener le produit en taille $n$ à 7 produits en taille $n/2$, plus un certain nombre C (ici $C = 18$) d'additions en taille $n/2$, utilisées pour fabriquer les combinaisons linéaires des blocs de A et de X, et celles des produits $q_i$.

La complexité $S(n)$ de l'algorithme de Strassen satisfait donc à la récurrence

$$S(n) \le 7S\left(\frac{n}{2}\right) + C\left(\frac{n}{2}\right)^2.$$

En invoquant le lemme « diviser pour régner » avec les choix $m = 7, p = s = 2, \kappa = 1$ et $q = 4$, on déduit l'inégalité

$$S(n) \le \left(1 + \frac{C}{3}\right) n^{\log_2 7},$$
qui prouve la partie (d) du Théorème 8.2.

L'idée de base de l'algorithme de Strassen est tout à fait générale : améliorer la multiplication en petite taille permet d'améliorer l'exposant.

**Exercice 8.2** Supposons que, pour un certain entier $r$ fixé, on sache multiplier deux matrices de taille $r$ à coefficients dans un corps $\mathbb{K}$ en $r^\theta$ multiplications non commutatives dans $\mathbb{K}$ et un nombre arbitraire d'additions dans $\mathbb{K}$. Montrer que $\theta$ est un exposant réalisable pour la multiplication des matrices à coefficients dans $\mathbb{K}$. ■

En pratique, il peut être intéressant de ne pas descendre récursivement jusqu'à la multiplication des matrices $1 \times 1$, mais d'arrêter la récursion à une taille $s > 1$ et d'utiliser l'algorithme naïf (ou tout autre algorithme, même commutatif !) en cette taille. Cela permet d'optimiser la constante du $O(\cdot)$.

**Exercice 8.3** Soit $n$ une puissance de 2. Établir un algorithme hybride de multiplication matricielle en taille $n \times n$, qui fait appel à l'algorithme de Strassen pour $n \ge 2^d$ et à l'algorithme naïf pour $n < 2^d$. Montrer que la complexité $MM(n)$ de cet algorithme vérifie $MM(n) \le \mu(d)n^{\log_2 7}$ pour tout $n \ge 2^d$, où $\mu(d)$ est une fonction qui dépend uniquement de $d$. Trouver la valeur de $d$ qui minimise $\mu(d)$. ■

Lorsque la taille $n$ des matrices à multiplier n'est pas une puissance de 2, on peut rajouter artificiellement des lignes et des colonnes aux deux matrices de façon à se ramener au cas où $n$ est une puissance de 2.

**Exercice 8.4** Écrire une version récursive de l'algorithme de Strassen valable en taille arbitraire $n$ (non nécessairement une puissance de 2), et qui procède en rajoutant une ligne et une colonne lorsque $n$ est impair. Montrer que la complexité de cet algorithme vérifie $MM(n) \le C(n)n^{\log_2 7}$ et comparer $C(n)$ à ce que l'on obtient en remplaçant $n$ par la puissance de 2 supérieure à $n$. ■

L'algorithme de Strassen s'applique aussi au cas où le corps de base $\mathbb{K}$ est remplacé par un anneau $\mathbb{A}$. Lorsque cet anneau $\mathbb{A}$ est une extension de $\mathbb{K}$, la question est d'estimer sa complexité en termes de nombre d'opérations dans $\mathbb{K}$.

**Exercice 8.5** Soit $\mathbb{K}$ un corps et soient $d$ et $n$ deux entiers naturels strictement positifs. Estimer la complexité de l'algorithme de multiplication de Strassen appliqué à deux matrices carrées de taille $n$, dont les éléments sont des polynômes de degré au plus $d$ de $\mathbb{K}[X]$. ■

Plus généralement, tout algorithme de multiplication dans $M_n(\mathbb{K})$ utilisant $MM(n) = O(n^\theta)$ opérations dans $\mathbb{K}$ fournit par application à $\mathbb{A} = \mathbb{K}[X]/(X^{2d+1})$ un algorithme de multiplication de matrices polynomiales de degré $d$ de complexité $MM(n,d) = O(n^\theta M(d))$ opérations de $\mathbb{K}$. Un algorithme de meilleure complexité permettant d'obtenir $MM(n,d) = O(n^\theta d + n^2 M(d))$ sera présenté au Chapitre 5.

### Interprétation des formules de Strassen

Il n'est pas immédiat de donner une intuition aux formules de Strassen (à la différence de celle de Karatsuba, qui repose de manière cachée sur l'évaluation de polynômes en 0, 1, +∞). Nous allons présenter un raisonnement qui permet de motiver
les formules de Strassen et de les généraliser. L'explication appartient à Fiduccia ; elle est postérieure à la découverte de l'algorithme<fnref n="3" />.

Le point de départ est l'observation que la recherche d'un algorithme non commutatif pour effectuer le produit de matrices $R = AX$ avec

$$A = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \quad \text{et} \quad X = \begin{pmatrix} x & y \\ z & t \end{pmatrix}$$

revient à donner un algorithme pour la multiplication matrice-vecteur $Mv$, où

$$M = \begin{pmatrix} a & b & 0 & 0 \\ c & d & 0 & 0 \\ 0 & 0 & a & b \\ 0 & 0 & c & d \end{pmatrix} \quad \text{et} \quad v = \begin{pmatrix} x \\ z \\ y \\ t \end{pmatrix}.$$

La suite de l'argument repose sur la remarque suivante : toute matrice (de taille arbitraire) de la forme

$$\begin{pmatrix} \alpha & \alpha \\ \varepsilon\alpha & \varepsilon\alpha \end{pmatrix}, \quad \begin{pmatrix} \alpha & -\alpha \\ \varepsilon\alpha & -\varepsilon\alpha \end{pmatrix} \quad \text{ou} \quad \begin{pmatrix} \alpha & \varepsilon\alpha \\ -\alpha & -\varepsilon\alpha \end{pmatrix},$$

où $\varepsilon \in \{0, 1\}$ et $\alpha \in \mathbb{K}$, dont tous les éléments, sauf les quatre marqués explicitement, sont nuls, peut être multipliée par un vecteur en une seule multiplication dans $\mathbb{K}$. On appelle une telle matrice *élémentaire*.

L'idée est alors de décomposer la matrice M en somme de 7 matrices élémentaires. On commence par soustraire à M les deux matrices élémentaires

$$E_1 = \begin{pmatrix} a & a \\ a & a \end{pmatrix} \quad \text{et} \quad E_2 = \begin{pmatrix} d & d \\ d & d \end{pmatrix},$$

correspondant aux éléments diagonaux de A. La matrice ainsi obtenue est

$$M_1 = \begin{pmatrix} b-a & & & \\ c-a & d-a & & \\ & a-d & b-d & \\ & c-d & & \end{pmatrix}.$$

<footnote n="3">3. Strassen lui-même ne se souvient pas exactement de l'origine de ses formules. Selon Landsberg, « Strassen was attempting to prove, by process of elimination, that such an algorithm did not exist when he arrived at it ». En tout cas, Strassen affirme « First I had realized that an estimate tensor rank < 8 for two by two matrix multiplication would give an asymptotically faster algorithm. Then I worked over $\mathbb{Z}/2\mathbb{Z}$ — as far as I remember — to simplify matters ». Autrement dit, il y a dû y avoir recherche « à la main ».</footnote>
Par construction, elle admet la propriété que l'élément (1,1) du deuxième bloc est l'opposé de l'élément (2,2) du premier. En soustrayant à $M_1$ la matrice élémentaire

$$E_3 = \begin{pmatrix} d-a & a-d \\ d-a & a-d \end{pmatrix}$$

on obtient la matrice

$$M_2 = \begin{pmatrix} c-a & b-a & & \\ & a-d & d-a & \\ & & c-d & b-d \end{pmatrix},$$

qui se décompose comme $M_3 + M_4$, avec

$$M_3 = \begin{pmatrix} b-a \\ a-d & b-d \end{pmatrix} \quad \text{et} \quad M_4 = \begin{pmatrix} c-a & d-a \\ & c-d \end{pmatrix}.$$

Or, puisque $a-d = (b-d) - (a-b)$ et $d-a = (c-a) - (d-b)$, chacune des matrices $M_3$ et $M_4$ s'écrit comme somme de 2 matrices élémentaires

$$M_3 = E_4 + E_5, \quad \text{où} \quad E_4 = \begin{pmatrix} b-a \\ a-b \end{pmatrix}, \quad E_5 = \begin{pmatrix} b-d & b-d \end{pmatrix},$$

et

$$M_4 = E_6 + E_7, \quad \text{où} \quad E_6 = \begin{pmatrix} c-a & c-a \end{pmatrix}, \quad E_7 = \begin{pmatrix} d-c \\ c-d \end{pmatrix}.$$

En revenant à l'égalité $^t(c r_{1,1} r_{2,1} r_{1,2} r_{2,2}) = (\sum_{i=1}^7 E_i)^t v$ traduisant $R = AX$, on obtient finalement les formules données dans les étapes (3) et (4) de l'Algorithmme 8.1.

**Exercice 8.6** Généraliser l'algorithme de Strassen aux matrices de taille arbitraire. Plus précisément, étendre la preuve des formules de Strassen pour le produit matriciel en taille $2 \times 2$ afin de montrer que le produit de deux matrices carrées de taille $n$ peut s'effectuer en $n^3 - n(n-1)/2$ multiplications scalaires non commutatives.

Est-ce que cela permet d'améliorer la borne de Strassen $\log_2 7$ sur l'exposant de l'algèbre linéaire ? ■
### Peut-on faire mieux que $2,81$ ?

Strassen ayant prouvé que l'exposant $\omega_{\text{mul}}$ du produit de matrices était sous-cubique, la question naturelle est devenue de savoir quelle est sa véritable valeur.

Pour ce faire, depuis les années 1970, de nombreux outils ont été développés, les uns plus sophistiqués que les autres, et l'on fit preuve de beaucoup d'ingéniosité. Des notions nouvelles sont apparues, comme celles de complexité bilinéaire, de rang tensoriel et d'algorithmes approchés, utilisées de manière intensive notamment par Bini, Pan, Schönhage, Strassen, Winograd et d'autres. Ces notions ne seront pas abordées dans cet ouvrage.

Le meilleur algorithme actuel, datant de 2014 et dû à François Le Gall, est de complexité $O(n^{2,3728639})$. Il est cependant conjecturé que $\omega_{\text{mul}} = 2$, mais pour l'instant la preuve (ou la réfutation) de cette conjecture semble hors de portée.

### En pratique

Les algorithmes rapides pour l'algèbre linéaire ont longtemps eu mauvaise presse; concrètement, l'algorithme de Strassen, et, de manière plus marginale, un algorithme d'exposant 2,77 dû à Pan et repris par Kaporin sont les seuls algorithmes « rapides » (avec un exposant inférieur à 3) qui ont été implantés avec succès.

Dans une bonne implantation, disons sur un corps fini défini modulo un entier de taille « raisonnable », les algorithmes de Winograd et de Waksman sont immédiate-ment rentables. L'algorithme de Strassen devient ensuite meilleur pour des tailles de l'ordre de quelques dizaines (64 est une borne raisonnable).

L'immense majorité des autres algorithmes connus reposent sur des techniques très complexes, qui se traduisent par la présence d'énormes constantes (et de facteurs logarithmiques) dans leurs estimations de complexité. En conséquence, dans leurs versions actuelles, il ne peuvent pas être rentables pour des tailles inférieures à des millions ou des milliards ...

## 8.3 Autres problèmes d'algèbre linéaire

Il est important de savoir multiplier les matrices, mais il est encore plus utile de les inverser, de calculer leur polynôme caractéristique, etc.

### Élimination de Gauss

Rappelons que l'algorithme classique d'élimination de Gauss permet de résoudre la plupart des problèmes classiques en algèbre linéaire de manière bien plus efficace que l'algorithme naïve.

**Définition 8.3 — Forme échelonnée.** Une matrice est échelonnée en ligne lorsque

1. les lignes nulles sont toutes en dessous des lignes non nulles,

2. le premier coefficient non nul (appelé *pivot*) de chaque ligne non nulle est strictement à droite du premier coefficient non nul de la ligne précédente.

Une forme échelonnée en ligne d'une matrice A est une matrice échelonnée en ligne B équivalente à gauche à A (c'est-à-dire telle qu'il existe une matrice carrée inversible P avec $A = PB$).
L’algorithme de Gauss (sans pivot sur les colonnes) calcule, en n’utilisant que des opérations élémentaires sur les lignes (qui correspondent à des multiplications à gauche par des matrices inversibles), une forme échelonnée en ligne de A, sur laquelle le rang et le déterminant de A se lisent directement.

En effet, les lignes non nulles d’une forme échelonnée B forment une base de l’espace vectoriel engendré par les lignes de la matrice de départ A. Le nombre de lignes non nulles de B est égal au rang de A. Si A est carrée, son déterminant est égal au produit des éléments diagonaux de B.

Une variante de l’algorithme de Gauss, dite de Gauss–Jordan, produit même une
forme échelonnée réduite de A : il s’agit d’une forme échelonnée de A, dont les pivots
valent 1, les autres coefficients dans les colonnes des pivots étant nuls. Appliquée à
la matrice augmentée $\tilde{A} = (A|I_n)$ obtenue par concaténation de la matrice A et de la
matrice identité $I_n$, cette variante est utilisée dans le calcul de l’inverse $A^{-1}$. En effet,
l’algorithme de Gauss–Jordan transforme la matrice $\tilde{A}$ en une matrice équivalente
dont le bloc gauche est l’identité, c’est-à-dire qu’il remplace $\tilde{A}$ par la matrice $(I_n|A^{-1})$.
Le même algorithme permet de résoudre le système $Ax = b$, où $b \in K^n$, en bordant la
matrice A par le vecteur $b$.

Une autre variante de l’algorithme d’élimination de Gauss permet de calculer
une décomposition LUP qui, à son tour, permet de résoudre des systèmes linéaires,
d’inverser des matrices, de calculer des déterminants, etc.

**Définition 8.4 — Matrices L, U, P et décomposition LUP.** Une matrice est dite *triangulaire supérieure*, resp. *triangulaire inférieure*, si les éléments situés en dessous (resp. au dessus) de la diagonale principale sont nuls. Une *matrice de permutation* est une matrice carrée P dont les coefficients valent 0 ou 1, et dont chaque ligne et colonne ne contient qu'un seul élément non nul. Une *décomposition LU*, resp. *LUP*, d'une matrice A est une factorisation de la forme $A = LU$, resp. $A = LUP$, avec L triangulaire inférieure, U triangulaire supérieure et P une matrice de permutation.

En résumé, nous allons admettre l’énoncé suivant.

**Théorème 8.3 — Élimination gaussienne.** Pour toute matrice $A \in M_n(\mathbb{K})$, il est
possible de calculer en $O(n^3)$ opérations dans $\mathbb{K}$ :

1. le rang $rg(A)$, le déterminant $det(A)$, et l'inverse $A^{-1}$ si A est inversible;

2. une base (affine) de solutions de $Ax = b$, pour tout $b$ dans $K^n$;

3. une décomposition LUP, et une forme échelonnée réduite de A.

Résultat principal

**Théorème 8.4 —** L’élimination gaussienne n’est pas optimale. Soit $\mathbb{K}$ un corps et soit $\theta$ un exposant réalisable pour la multiplication des matrices à coefficients dans $\mathbb{K}$. Pour toute matrice $A$ de $M_n(\mathbb{K})$, il est possible de calculer :

a. le déterminant $det(A)$;

b. l'inverse $A^{-1}$ (si A est inversible);

c. la solution $x \in K^n$ du système $Ax = b$, pour tout $b \in K^n$, si A est inversible;
d. le polynôme caractéristique de A;

e. une décomposition LUP de A;

f. le rang $rg(A)$ et une forme échelonnée réduite de A;

g. une base du noyau de A;

en $\tilde{O}(n^\theta)$ opérations dans $\mathbb{K}$.

La preuve des parties (e), (f) et (g) dépasse le cadre de cet ouvrage. En termes d'expô-
sants, les assertions (a) et (d) impliquent les inégalités $\omega_{inv} \le \omega_{mul}$ et $\omega_{polcar} \le \omega_{mul}$.
L'inégalité $\omega_{det} \le \omega_{polcar}$ est évidente, puisque le terme constant du polynôme carac-
téristique de A est égal à $(-1)^n \det(A)$.

Dans la suite, nous allons prouver (a)–(d), sous des hypothèses simplificatrices. Nous montrerons également que $\omega_{mul} \le \omega_{inv}$ et laisserons en exercice (Exercice 8.13) la preuve de $\omega_{inv} \le \omega_{det}$. Utilisées conjointement, ces inégalités permettent de prouver le Théorème 8.1.

La multiplication n'est pas plus difficile que l'inversion

Soient $A$ et $B$ des matrices $n \times n$. On souhaite calculer $C = AB$. Pour cela, on pose

$$
D = \begin{pmatrix} I_n & A & 0 \\ 0 & I_n & B \\ 0 & 0 & I_n \end{pmatrix}.
$$

On a alors l'identité remarquable

$$
D^{-1} = \begin{pmatrix} I_n & -A & AB \\ 0 & I_n & -B \\ 0 & 0 & I_n \end{pmatrix},
$$

qui permet donc de ramener le produit en taille *n* à l'inversion en taille *3n*. Cela prouve notre affirmation et l'inégalité

$$
\omega_{\text{mul}} \leq \omega_{\text{inv}}.
$$

Exercice 8.7 Soit $\mathbb{K}$ un corps et soit $T(n)$ la complexité du produit de deux matrices triangulaires inférieures de taille $n \times n$ et à coefficients dans $\mathbb{K}$. Montrer qu'il est possible de multiplier deux matrices arbitraires de $M_n(\mathbb{K})$ en $O(T(n))$ opérations dans $\mathbb{K}$. ■

L'inversion, le calcul de déterminant et la résolution de système ne sont pas
plus difficiles que la multiplication.

Nous allons montrer qu'il est possible d'inverser des matrices $n \times n$ par un algorithme de type « diviser pour régner » en $\tilde{O}(n^\theta)$ opérations arithmétiques, quel que soit l'exposant réalisable $\theta$.

Inverse

L'algorithme d'inversion matricielle présenté ici, dû à Strassen, est un algorithme récursif, qui peut s'interpréter comme un « pivot de Gauss par blocs ». Cet algorithme nécessite d'inverser certaines sous-matrices de A; afin de garantir sa correction, nous
**Entrée** Une matrice « générique » $A \in M_n(\mathbb{K})$, avec $n = 2^k$.

**Sortie** Son inverse $A^{-1}$.

1. Si $n = 1$, renvoyer $A^{-1}$.

2. Décomposer $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$, avec $a, b, c, d$ dans $M_{n/2}(\mathbb{K})$.

3. Calculer $e := a^{-1}$ récursivement.

4. Calculer $Z := d - ceb$.

5. Calculer $t := Z^{-1}$ récursivement.

6. Calculer $y := -ebt$, $z := -tce$, $x := e + ebtce$.

7. Renvoyer $\begin{pmatrix} x & y \\ z & t \end{pmatrix}$.

**Algorithme 8.2** – Algorithme de Strassen pour inverser une matrice.

allons faire l'hypothèse simplificatrice que *toutes ces matrices sont inversibles*. Le traitement du cas général (A matrice arbitraire) est plus délicat, et passe par le calcul efficace, toujours de type « diviser pour régner », d'une décomposition LUP de A.

Le point de départ est l'identité suivante (non commutative!), dans laquelle on suppose $n = 2$ et $a, Z \in \mathbb{K} \setminus \{0\}$ :

$$A = \begin{pmatrix} a & b \\ c & d \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ ca^{-1} & 1 \end{pmatrix} \times \begin{pmatrix} a & 0 \\ 0 & Z \end{pmatrix} \times \begin{pmatrix} 1 & a^{-1}b \\ 0 & 1 \end{pmatrix}, \quad (8.2)$$

et où $Z = d - ca^{-1}b$ est le complément de Schur de $a$ dans A.

Cette identité se déduit aisément grâce à la méthode du pivot de Gauss appliquée à A et permet d'obtenir la factorisation suivante

$$ \begin{align} \begin{pmatrix} a & b \\ c & d \end{pmatrix}^{-1} &= \begin{pmatrix} 1 & -a^{-1}b \\ 0 & 1 \end{pmatrix} \times \begin{pmatrix} a^{-1} & 0 \\ 0 & Z^{-1} \end{pmatrix} \times \begin{pmatrix} 1 & 0 \\ -ca^{-1} & 1 \end{pmatrix} \tag{8.3} \\ &= \begin{pmatrix} a^{-1} + a^{-1}bZ^{-1}ca^{-1} & -a^{-1}bZ^{-1} \\ -Z^{-1}ca^{-1} & Z^{-1} \end{pmatrix}. \notag \end{align} $$

Le point important est que l'identité (8.3) est « non commutative » : elle reste valable dans tout anneau non commutatif, pourvu que les éléments $a$ et $Z = d - ca^{-1}b$ soient inversibles. Elle permet donc une application récursive. L'algorithme 8.2 résume la méthode d'inversion. Dans chaque appel récursif l'algorithme effectue deux inversions, ainsi qu'un certain nombre de produits. Plus précisemment, avec les notations de l'Algorithme 8.2, les 6 produits $ce, bt, ceb = (ce)b, ebt = e(bt), (ebt)(ce)$ et $t(ce)$ sont suffisants. Cette analyse mène au résultat suivant.

**Théorème 8.5** Si toutes les sous-matrices rencontrées au cours de l'Algorithme 8.2 sont inversibles, le coût de l'inversion de A est de $3MM(n) + O(n^2)$.

**Démonstration.** L'algorithme ci-dessus montre que la complexité $l(n)$ de l'inversion
satisfait à la récurrence

$$l(n) \leq 2l\left(\frac{n}{2}\right) + 6MM\left(\frac{n}{2}\right) + Cn^2,$$

C étant une constante. Il suffit alors d'appliquer le lemme « diviser pour régner » avec les choix $m = p = s = 2$, $\kappa = 1$ et $q = 4$. ■

Une première conséquence immédiate du Théorème 8.5 est que la résolution du système linéaire $Ax = b$, pour $b \in \mathbb{K}^n$, et A inversible vérifiant les hypothèses du Théorème 8.5, peut s'effectuer également en $O(\text{MM}(n)) = O(n^\theta)$ opérations dans $\mathbb{K}$.

## Déterminant

La factorisation (8.2) montre qu'une légère adaptation de l'algorithme d'inversion permet de calculer le déterminant de A en même temps que son inverse, en la même complexité. Pour cela, il suffit de remplacer les étapes (1), (3), (5) et (7) de l'Algorithme 8.2 par les étapes (1'), (3'), (5') et (7') décrites ci-dessous.

1'. Si $n = 1$, renvoyer $A^{-1}$ et A.

3'. Calculer $e := a^{-1}$ et $d_a := \det(a)$ récursivement.

5'. Calculer $t := Z^{-1}$ et $d_Z := \det(Z)$ récursivement.

7'. Renvoyer $\begin{pmatrix} x & y \\ z & t \end{pmatrix}$ et $d_a d_Z$.

## Affaiblissement des hypothèses

À ce stade, nous avons démontré les parties (a)–(c) du Théorème 8.4 sous une hypothèse supplémentaire. Dans le cas où $\mathbb{K}$ est un sous-corps du corps des réels $\mathbb{R}$, on peut s'affranchir de cette hypothèse, en observant que, quelle que soit la matrice A inversible dans $M_n(\mathbb{R})$, la matrice $B = {}^tA \cdot A$ est définie positive et vérifie les hypothèses du Théorème 8.5, et donc l'inverse de A peut se calculer grâce à l'égalité $A^{-1} = B^{-1} \cdot {}^tA$ en $O(\text{MM}(n)) = O(n^\theta)$ opérations dans $\mathbb{K}$.

Il n'est pas difficile de prouver que si la matrice A est de type L (triangulaire inférieure) ou U (triangulaire supérieure) alors l'algorithme de Strassen calcule son inverse sans avoir besoin d'aucune hypothèse supplémentaire. Cette remarque est à la base d'un algorithme général d'inversion dû à Bunch et Hopcroft. Cet algorithme calcule d'abord une décomposition LUP d'une matrice inversible arbitraire<fnref n="4" /> *sur un corps quelconque* $\mathbb{K}$ en $O(\text{MM}(n))$ opérations, et arrive à en déduire le calcul d'inverse (et aussi la résolution de système) pour le même prix.

L'algorithme d'inversion de Strassen est en outre à la base de toute une algorithmique pour les matrices structurées, qui sera présentée au Chapitre 10.

## Résolution des sytèmes surdéterminés

Lorsqu'une matrice n'est pas carrée, le calcul de son noyau se ramène aisément au cas carré, comme illustré dans le résultat suivant qui nous sera utile dans le Chapitre 22 sur la factorisation des polynômes à plusieurs variables.

<footnote n="4">4. La détection de l'inversibilité est également réalisable en $O(\text{MM}(n))$ opérations, grâce à un algorithme qui calcule le rang de A en cette complexité; cet algorithme dépasse notre cadre.</footnote>
**Théorème 8.6** Soit $\mathbb{K}$ un corps et supposons qu'on dispose d'un algorithme de multiplication dans $M_n(\mathbb{K})$ de complexité $\text{MM}(n) = O(n^\theta)$, avec $\theta \ge 2$. Le calcul d'une base du noyau d'une matrice $A$ à $n$ colonnes et $m \ge n$ lignes sur $\mathbb{K}$ coûte $\tilde{O}(mn^{\theta-1})$ opérations dans $\mathbb{K}$.

*Démonstration.* Le Théorème 8.4 implique ce lemme si $m = n$. Donc, quitte à ajouter des lignes nulles à la matrice $A$, on peut supposer que $m$ est un multiple de $n$. Notons $A_1$ la matrice composée des $m$ premières lignes de $A$, $A_2$ celle des $m$ suivantes, etc. On peut calculer une base du noyau de $A_1$ avec un coût $\tilde{O}(n^\theta)$ : notons $N_1$ la matrice dont les colonnes sont les vecteurs de cette base, de sorte que $A_1 N_1 = 0$. On calcule ensuite la matrice $N_2$ dont les colonnes forment une base du noyau de $A_2 N_1$, de sorte que $A_2 N_1 N_2 = 0$. On continue de la même manière avec $A_3$, en calculant $N_3$ telle que $A_3 N_1 N_2 N_3 = 0$, etc. Finalement les colonnes de $N_1 \cdots N_{m/n}$ forment une base du noyau de $A$. On obtient ainsi le coût total annoncé $\tilde{O}(mn^{\theta-1})$. ■

## Calcul du polynôme caractéristique

Soit $A$ une matrice de $M_n(\mathbb{K})$. Dans cette section, nous présentons un algorithme efficace, dû à Keller-Gehrig, pour le calcul du polynôme caractéristique $\chi_A(X) = \det(XI_n - A)$ de $A$. Nous ne décrirons en détail qu'un cas particulier très simple (et très fréquent). Plus exactement, nous ferons dans la suite l'hypothèse que $\chi_A(X) = X^n + p_{n-1}X^{n-1} + \cdots + p_0$ est irréductible dans $\mathbb{K}[X]$; en particulier, il coïncide avec le polynôme minimal de $A$. L'algorithme repose de façon cruciale sur la propriété suivante.

**Lemme 8.7** Soit $A$ une matrice de $M_n(\mathbb{K})$, dont le polynôme caractéristique est irréductible dans $\mathbb{K}[X]$. Soit $v$ un vecteur de $\mathbb{K}^n \setminus \{0\}$ et soit $P \in M_n(\mathbb{K})$ la matrice dont la $j$-ième colonne vaut $A^{j-1}v$, pour $1 \le j \le n$. Alors $P$ est inversible et la matrice $P^{-1}AP$ est de type compagnon.

Rappelons qu'une matrice $C = (c_{i,j})_{i,j=1}^n$ de $M_n(\mathbb{K})$ est appelée *de type compagnon* si ses $n-1$ premières colonnes ne contiennent que des éléments nuls, à l'exception des éléments $c_{2,1}, c_{3,2}, \dots, c_{n,n-1}$ qui valent tous 1. L'intérêt du Lemme 8.7 vient du fait classique que le polynôme caractéristique de la matrice $C$ vaut $X^n - \sum_{i=1}^n c_{i,n} X^{i-1}$ et n'a donc besoin d'aucune opération arithmétique pour être calculé.

*Démonstration.* Pour prouver ce lemme, on commence par observer que, grâce à l'hypothèse sur $A$, la famille $B = \{v_1, v_2, \dots, v_n\}$, où $v_i = A^{i-1}v$, forme une base du $\mathbb{K}$-espace vectoriel $\mathbb{K}^n$, et donc $P$ est inversible. En effet, en supposant le contraire, on déduit l'existence d'un polynôme non nul $Q \in \mathbb{K}[X]$ de degré strictement inférieur à $n$, tel que $Q(A)v = 0$. Sans perte de généralité, ce polynôme peut être supposé de degré minimal avec cette propriété. Comme par ailleurs $\chi_A(A)v = 0$, il s'ensuit que $Q$ divise $\chi_A$; or, ce dernier étant supposé irréductible, cela entraîne que $Q$ est constant, et donc $v = 0$, une contradiction.

Puisque $Av_{i-1} = v_i$, la matrice $C$ de l'application linéaire $w \mapsto Aw$ dans la base $B$ est de type compagnon. Le théorème du changement de base fournit alors l'égalité $C = P^{-1}AP$ et permet donc de conclure. ■
**Entrée** Une matrice A de $M_n(\mathbb{K})$, avec $n = 2^k$.

**Sortie** Son polynôme caractéristique $\det(XI_n - A)$.

1. Choisir $v$ dans $\mathbb{K}^n \setminus \{0\}$.

2. $M := A$; $P := v$.

3. Pour $i$ de 1 à $k$, remplacer $P$ par la concaténation horizontale de $P$ et $MP$, puis $M$ par $M^2$;

4. $C := P^{-1}AP$.

5. Renvoyer $(X^n + p_{n-1}X^{n-1} + \cdots + p_0)$, où $^t(-p_0, \ldots, -p_{n-1})$ est la dernière colonne de C.

**Algorithme 8.3** – Calcul du polynôme caractéristique par l’algorithme de Keller-Gehrig.

Les matrices $A$ et $C = P^{-1}AP$ étant semblables, elles ont le même polynôme caractéristique. L'idée de l'algorithme de Keller-Gehrig est alors naturelle : construire la matrice $P$ du Lemme 8.7, pour ensuite déterminer la matrice $C$ à l'aide d'une multiplication et d'une inversion, et lire les coefficients du polynôme caractéristique de $A$ sur les éléments de la dernière colonne de $C$.

Du point de vue de la complexité, la seule étape coûteuse est la construction de la matrice $P$. En effet, la méthode directe qui consiste à calculer la suite (dite de Krylov) $v, Av, \dots, A^{n-1}v$ par multiplications successives d'un vecteur par la matrice $A$ a un coût cubique en $n$. La remarque cruciale est que l'on peut regrouper les produits matrice-vecteur en plusieurs produits matrice-matrice, de la façon suivante : on calcule en $O(n^2)$ opérations les vecteurs $v$ et $Av$, on détermine ensuite par exponentiation binaire les $O(\log n)$ matrices $A, A^2, A^4, A^8, \dots$ et on termine le calcul des colonnes de $P$ par le calcul des produits $A^2 \times (v | Av) = (A^2v | A^3v)$, puis $A^4 \times (v | \dots | A^3v) = (A^4v | \dots | A^7v)$, etc. Chaque produit de ce type est effectué à l'aide d'un produit matriciel en taille $n \times n$, en rajoutant artificiellement des colonnes nulles aux facteurs droits. Cette méthode conduit à l'Algorithmme 8.3. Nous venons de prouver la version suivante plus faible du point (d) du Théorème 8.4 :

**Théorème 8.8** Soit $\mathbb{K}$ un corps et supposons qu'on dispose d'un algorithme de multiplication dans $M_n(\mathbb{K})$ de complexité $\text{MM}(n) = O(n^\theta)$, avec $\theta \ge 2$. Si le polynôme caractéristique d'une matrice $A$ de $M_n(\mathbb{K})$ est irréductible dans $\mathbb{K}[X]$, l'Algorithmme 8.3 le calcule en $O(\text{MM}(n) \log n)$ opérations dans $\mathbb{K}$.

Le calcul de la suite de Krylov sera utilisé au Chapitre 9 dans l'algorithme de Wiedemann pour la résolution des systèmes creux, et de façon moins transparente, au Chapitre 11 dans l'algorithme de Storjohann pour la résolution de systèmes linéaires à coefficients polynomiaux.

## Algorithmes sans division

Il est parfois utile de calculer le polynôme caractéristique sans effectuer de division dans l'anneau des coefficients de la matrice. Nous utiliserons en Partie V le résultat suivant, que nous admettrons ici : en comparaison avec le Théorème 8.4, l'évitement
des divisions se fait en perdant 1 sur l'exposant de complexité.

**Théorème 8.9** — **Berkowitz.** Le polynôme caractéristique d'une matrice carrée de taille $n$ à coefficients dans un anneau commutatif $\mathbb{A}$ peut se calculer en $O(n^{\theta+1})$ opérations dans $\mathbb{A}$, où $\theta$ est un exposant réalisable pour le produit.

Il est cependant aisé d'obtenir cette même borne de complexité en supposant les divisions par $2, \dots, n$ possibles dans l'anneau des coefficients (ce qui est le cas en présence d'une $\mathbb{Q}$-algèbre par exemple). En effet, la $i$-ième somme de Newton (définie en Section 3.3) du polynôme caractéristique de la matrice $\mathbb{A}$ s'obtient comme $\text{Tr}(\mathbb{A}^i)$. Il suffit donc d'utiliser l'Algorithmhe 3.5 pour obtenir les coefficients du polynôme caractéristique de $\mathbb{A}$ à partir de ces traces pour $i = 0, \dots, n-1$. Il s'agit de la *méthode de Le Verrier*.

## Exercices

**Exercice 8.8** — **Multiplication de matrices rectangulaires.** Si $\mathbb{A}$ et $\mathbb{B}$ sont deux matrices $n \times r$ et $r \times k$ avec $n \ge r$ et $k \ge r$, à coefficients dans un corps $\mathbb{K}$, expliquer comment calculer le produit $\mathbb{A}\mathbb{B}$ en $O(nkr^{\theta-2})$ opérations arithmétiques dans $\mathbb{K}$. ■

**Exercice 8.9** — **Multiplication et division rapides de nombres complexes.**

1. Montrer qu'il est possible d'effectuer le produit matrice-vecteur

$$ \begin{pmatrix} a & -b \\ b & a \end{pmatrix} \times \begin{pmatrix} x \\ y \end{pmatrix} $$

en utilisant seulement 3 multiplications scalaires.

2. En déduire que deux nombres complexes peuvent être multipliés en utilisant 3 multiplications réelles.

3. Montrer qu'on peut calculer le quotient de deux nombres complexes en utilisant au plus 7 multiplications et divisions réelles.

4. Donner un algorithme qui utilise seulement 6 opérations de ce type. ■

**Exercice 8.10** — **Multiplication rapide de quaternions.**

1. Soit $\mathbb{A} = (a_{ij})$ une matrice $n \times n$ à coefficients dans $\mathbb{K}$, telle que $a_{ij} = \pm a_{ji}$ pour $i < j$. Soit $v$ un vecteur quelconque de $\mathbb{K}^n$. Si $\ell$ est le nombre d'éléments non nuls de l'ensemble $\{a_{ij}, i < j\}$, montrer que $n + \ell$ multiplications suffisent pour calculer le produit $\mathbb{A}v$.

2. Montrer qu'on peut calculer le produit d'une matrice $n \times n$ et d'une matrice $n \times 2$ en utilisant $\frac{3}{2}n^2 + \frac{5}{2}n$ multiplications scalaires. Indication : Remplacer le produit $\mathbb{M} \times (v_1 | v_2) \text{ par } \begin{pmatrix} 0 & \mathbb{M} \\ \mathbb{M} & 0 \end{pmatrix} \times \begin{pmatrix} v_1 \\ v_2 \end{pmatrix}$, puis décomposer la matrice $\mathbb{A} = \begin{pmatrix} 0 & \mathbb{M} \\ \mathbb{M} & 0 \end{pmatrix}$ en somme de deux matrices vérifiant les hypothèses de la question précédente.

3. Un *quaternion* est une combinaison formelle $a \cdot 1 + b \cdot i + c \cdot j + d \cdot k$, avec $a, b, c, d \in \mathbb{R}$
et $i, j, k$ éléments vérifiant la table de multiplication interne

<table><thead><tr><th>·</th><th>1</th><th>i</th><th>j</th><th>k</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>i</td><td>j</td><td>k</td></tr><tr><th>i</th><td>i</td><td>-1</td><td>k</td><td>-j</td></tr><tr><th>j</th><td>j</td><td>-k</td><td>-1</td><td>i</td></tr><tr><th>k</th><td>k</td><td>j</td><td>-i</td><td>-1</td></tr></tbody></table>

(par exemple, $i \cdot j = k$ et $k \cdot j = -i$). La multiplication $\cdot$ se prolonge par $\mathbb{R}$-linéarité à l'ensemble $\mathbb{H}$ des quaternions et en fait un corps non commutatif.

Montrer que deux quaternions arbitraires peuvent être multipliés en utilisant seulement 10 multiplications réelles.

■

**Exercice 8.11** Soit P un polynôme de degré au plus $n$ à coefficients dans un corps $\mathbb{K}$. Soit $\theta > 2$ un réel tel que $O(n^{\theta})$ opérations dans $\mathbb{K}$ suffisent pour multiplier deux matrices arbitraires de $M_n(\mathbb{K})$.

1. Donner un algorithme pour évaluer $P$ sur $\lceil \sqrt{n} \rceil$ éléments de $\mathbb{K}$ en $O(n^{\theta/2})$ opérations.

2. Si $A \in M_n(\mathbb{K})$, montrer qu'on peut calculer la matrice $P(A)$ en $O(n^{\theta+1/2})$ opérations de $\mathbb{K}$.

3. Si $Q \in \mathbb{K}[X]$ est un autre polynôme de degré au plus $n$, montrer qu'il est possible de calculer les $n$ premiers coefficients du polynôme $P(Q(X))$ en $O(n^{\theta+1/2})$ opérations dans $\mathbb{K}$.

*Indication :* Écrire $P$ sous la forme $P_0(X) + P_1(X)X^d + P_2(X)(X^d)^2 + \dots$, avec $d$ bien choisi et $P_i(X)$ de degrés au plus $d-1$.

■

**Exercice 8.12** — Inversion de Strassen des matrices polynomiales. Soit $M(X)$ une matrice de $M_n(\mathbb{K}[X])$ une matrice polynomiale carrée de taille $n$ et de degré borné par $d$. Exhiber une borne supérieure pour la complexité (exprimée en termes d'opérations arithmétiques dans $\mathbb{K}$ et en fonction des deux paramètres $n$ et $d$) de l'algorithme de Strassen pour le calcul de $M(X)^{-1}$, sous l'hypothèse que toutes les matrices rencontrées au cours de l'algorithme sont inversibles.

■

**Exercice 8.13** — Dérivation automatique et calcul d’inverse. Le but de cet exercice est de démontrer que le calcul de l’inverse d’une matrice n’est pas plus difficile que le calcul du déterminant, en utilisant des techniques de dérivation automatique.

Pour cela, formalisons brièvement ce qu'on entend par *programme*. Dans notre contexte, un programme prend en entrée des indéterminées $Y_1, \dots, Y_m$ sur un corps $\mathbb{K}$, et effectue une suite d'instructions $g_1, \dots, g_L$. Pour tout $i = 1, \dots, L$, l'instruction $g_i$ consiste simplement en la définition d'un polynôme $G_i$, de la forme

$$g_i : G_i = \mathrm{Arg}_i \circ p_i \mathrm{Arg}'_i,$$

selon les règles suivantes :

— Arg$_i$ et Arg'$_i$ sont dans $\mathbb{K} \cup \{G_1, \dots, G_{i-1}\} \cup \{Y_1, \dots, Y_m\}$.
— op<sub>i</sub> est dans {+, -, x}.

La *taille* d'un tel programme est L. On dit qu'un tel programme *calcule* un polynôme F si F appartient à l'ensemble {G<sub>1</sub>, ..., G<sub>L</sub>} ; on dit de même que le programme calcule {F<sub>1</sub>, ..., F<sub>s</sub>} si tous les F<sub>i</sub> sont dans {G<sub>1</sub>, ..., G<sub>L</sub>}.

Par exemple, le programme ayant pour entrées Y₁, Y₂, Y₃, Y₄, pour instructions

$$
G_1 = Y_1 \times Y_4
$$

$$
G_2 = Y_2 \times Y_3
$$

$$
G_3 = G_1 - G_2
$$

calcule le déterminant Y₁Y₄ – Y₂Y₃ de la matrice

$$
\begin{pmatrix}
Y_1 & Y_2 \\
Y_3 & Y_4
\end{pmatrix};
$$

il a taille 3.

Si G est dans $\mathbb{K}[Y_1, \dots, Y_m]$, on appelle *gradient* de G l'ensemble de ses $m$ dérivées partielles $\partial G / \partial Y_i$, $i = 1, \dots, m$.

1. Soient G et H dans $\mathbb{K}[Y_1, \dots, Y_m]$. Exprimer les dérivées partielles de G+H, G-H et G×H, en fonction de G, H et des dérivées partielles de G et H.
En déduire qu'étant donné un programme de taille L, qui calcule un polynôme G de $\mathbb{K}[Y_1, \dots, Y_m]$, on peut en déduire un programme de taille $O(Lm)$ qui calcule G et son gradient.

2. Soient G dans $\mathbb{K}[Y_1, \dots, Y_m]$ et H dans $\mathbb{K}[Y_1, \dots, Y_{m+1}]$, tels que

$$
G = H(Y_1, \dots, Y_m, Y_i \times Y_j),
$$

avec 1 ≤ i, j ≤ m. Exprimer les dérivées partielles de G en fonction de H et de
ses dérivées partielles. Traiter ensuite les cas

$$
G = H(Y_1, \dots, Y_m, Y_i \pm Y_j).
$$

3. Soient G dans $\mathbb{K}[Y_1, \dots, Y_m]$ et H dans $\mathbb{K}[Y_1, \dots, Y_{m+1}]$, tels que

$$
G = H(Y_1, \dots, Y_m, Y_i \text{ op } Y_j),
$$

avec 1 ≤ i, j ≤ m, et op dans {+, -, x}. On suppose connu un programme de taille L qui calcule H et son gradient. Montrer qu'on peut en déduire un programme de taille L + O(1) qui calcule G et son gradient.

En déduire par récurrence que si un polynôme G de $\mathbb{K}[Y_1, \dots, Y_m]$ peut être calculé par un programme de taille L, G et son gradient peuvent être calculés par un programme de taille O(L), où la constante dans le O(·) ne dépend pas du nombre de variables m.
4. Soient $X_{1,1}, \dots, X_{1,n}, \dots, X_{n,1}, \dots, X_{n,n}$ des indéterminées, et soit D le déterminant de la matrice

$$ M = \begin{pmatrix} X_{1,1} & \cdots & X_{1,n} \\ \vdots & & \vdots \\ X_{n,1} & \cdots & X_{n,n} \end{pmatrix}. $$

Montrer que le gradient de D permet de retrouver les entrées de la comatrice N de M, où N est l'unique matrice satisfaisant

$$ M \cdot {}^t N = {}^t N \cdot M = D I_n, $$

$I_n$ étant la matrice identité de taille $n$, et ${}^t N$ la transposée de N.

5. Montrer qu'étant donné un programme de taille L qui calcule D, on peut en déduire un programme de taille O(L), effectuant éventuellement des divisions, qui calcule les coefficients de l'inverse de M.

6. En déduire que $\omega_{inv} \le \omega_{det}$.

■

## Notes

L'optimalité de l'algorithme naïf pour le produit matrice-vecteur mentionnée au début de la Section 8.2 a été prouvée par Winograd [Win67]. Le résultat de l'Exercice 8.9 est dû à Smith [Smi62]; Alt et van Leeuwen [AL81] et Lickteig [Lic87] ont prouvé que la borne 6 est optimale. Le résultat de l'Exercice 8.10 est dû à Fiducia [Fid72c]; la borne 10 peut être améliorée à 8, et celle-ci est optimale [BCS97, (17.32)].

Les parties (b) et (c) du Théorème 8.2 sont dues à Winograd [Win68] et Waksman [Wak70]. Paterson [Pat75] a prouvé que, si l'on s'interdit les soustractions, la multiplication matricielle cubique est optimale. Le Théorème 8.4, dont le titre est emprunté à l'article de Strassen de 1969 [Str69], montre que du point de vue de la complexité on peut faire mieux que la méthode du pivot de Gauss. C'est un résultat surprenant car Klyuyev et Kokovkin-Shcherbak [KK65] avaient prouvé que l'élimination gaussienne est optimale si l'on se restreint à des opérations sur des lignes et des colonnes. C'est surtout un résultat qui a ouvert la voie à tout un domaine de recherche toujours très actif, la théorie de la complexité bilinéaire. Les livres de Pan [Pan84b], de Bürgisser, Clausen et Shokrollahi [BCS97, §15–16] et de Abdeljaoued et Lombardi [AL04, §7–10], ainsi que les articles de synthèse de Pan [Pan84a] et de Strassen [Str90] constituent de bonnes références sur ce sujet.

L'interprétation des formules de Strassen donnée en Section 8.2 et l'Exercice 8.6 sont dus à Fiduccia [Fid72c]. L'algorithme d'inversion des matrices génériques en Section 8.3 est dû à Strassen [Str69]. La remarque de Schönhage permettant de traiter l'inversion rapide d'une matrice réelle arbitraire est tirée de [Sch73]. L'algorithme général d'inversion évoqué à la fin de la Section 8.3, ainsi que la preuve des égalités $\omega_{inv} = \omega_{LUP} = \omega_{mul}$ ont été donnés par Bunch et Hopcroft [BH74]. La preuve de l'inégalité $\omega_{inv} \le \omega_{det}$ esquissée dans l'Exercice 8.13 est due à Baur et Strassen [BS83]. La réduction présentée en Section 8.3 provient de l'article de Munro [Mun73]. La partie (f) du Théorème 8.4 est due à Keller-Gehrig [Kel85], et la partie (g) est due à
Bürgisser, Karpinski et Lickteig [BKL91]. Le fait que l'exposant $\omega_{\text{mul}}$ ne dépend que de la caractéristique du corps de base a été prouvé par Schönhage [Sch81].

L'algorithme décrit dans la Section 8.3 a été conçu par Keller-Gehrig [Kel85]; il s'agit d'une version accélérée d'un algorithme classique de Krylov et Danilevskii [Dan37; Kry31]. Il est possible d'adapter cet algorithme au cas où A est une matrice arbitraire. Plus précisément, l'algorithme de Keller-Gehrig calcule dans ce cas une base de $K^n$ dans laquelle la matrice de l'endomorphisme de multiplication par A est triangulaire par blocs, avec des blocs diagonaux de type compagnon. Autrement dit, il procède à une mise en forme de Frobenius faible. Un raffinement de cet algorithme, dû à Giesbrecht [Gie95], permet de calculer la forme de Frobenius forte (matrice diagonale par blocs de type compagnon) également en complexité $\tilde{O}(\text{MM}(n))$; cet algorithme permet aussi le calcul, en la même complexité, du polynôme minimal et de tous les autres facteurs invariants de A. Les algorithmes de Giesbrecht [Gie95] sont probabilistes; des versions déterministes ont été données par Storjohann [Sto01]. Keller-Gehrig [Kel85] propose également un algorithme différent de celui décrit dans la Section 8.3, calculant le polynôme caractéristique d'une matrice générique en seulement $O(\text{MM}(n))$ opérations; voir l'Exercice 33.11 en page 599. Une version (probabiliste) de même complexité et sans hypothèse de généricité a été proposée par Pernet et Storjohann [PS07].

Giesbrecht [Gie95] et Storjohann [Sto01] montrent que $\omega_{\text{Frobenius}} = \omega_{\text{mul}}$ et que l'évaluation d'un polynôme de degré $d$ en une matrice de $M_n(K)$ peut se faire en $\tilde{O}(\text{MM}(n) + d)$ opérations dans $K$, ce qui améliore le résultat (2) de l'Exercice 8.11. Cet exercice est inspiré des articles de Borodin et Munro [BM71] pour (1), de Paterson et Stockmeyer [PS73] pour (2) et de Brent et Kung [BK78] pour (3); le résultat de (1) sera amélioré au Chapitre 5 et celui de (3) au Chapitre 3.

À l'heure actuelle, on ne sait pas si les exposants des problèmes suivants

— résoudre un système;

— calculer le rang;

— décider si une matrice est inversible;

sont eux aussi égaux à $\omega_{\text{mul}}$. Grâce au Théorème 8.4, on sait que ces problèmes ne sont pas plus difficiles que la multiplication, mais peut-être sont-ils plus faciles...

Un autre problème ouvert est le calcul du polynôme caractéristique par un algorithme déterministe de complexité $O(\text{MM}(n))$.

L'algorithme de Berkowitz a été décrit pour la première fois en 1984 [Ber84].
L'algorithme de Le Verrier se trouve dans son article de 1840 [Le 40].

### Complexité des algorithmes commutatifs

Il est naturel de se demander s'il est possible d'obtenir des algorithmes à la Winograd plus efficaces, en regroupant les termes de (8.1) en paquets plus gros plutôt que deux par deux. La réponse est négative. Harter [Har72] a montré qu'il n'existe pas de généralisation de (8.1) de la forme

$$ (\ell|c) = \lambda(a_1, x_1, a_2, x_2, \dots) + \mu(a_1, x_1, a_2, x_2, \dots) + f(a_1, a_2, \dots) + g(x_1, x_2, \dots), $$

où $\lambda$ et $\mu$ sont des formes linéaires et $f$ et $g$ des fonctions quelconques. Ja'Ja [JaJ79] a montré plus généralement qu'aucun algorithme commutatif ne peut apporter un gain d'un facteur supérieur à 2 sur l'algorithme naïf, et en particulier ne peut pas
amideorer l'exposant 3. De ce point de vue, les algorithmes de Winograd et Waksman sont optimaux dans la classe des algorithmes commutatifs.

Algorithme de Strassen

Winograd [FP74] a donné une version de l'algorithme de Strassen utilisant seule-
ment $C = 15$ additions et soustractions (au lieu de $C = 18$). Cela permet de descendre
la constante dans le $O(n^{\log_2 7})$ de 7 à 6 pour l'Algorithme 8.1. Cette constante a été
abaissée à 3,92 par Fischer [Fis74].

Dans le cas où la taille *n* de la matrice n'est pas une puissance de 2, il vaut
mieux éviter d'appliquer l'algorithme de Strassen à la matrice augmentée en taille la
puissance de 2 supérieure à *n*, comme le montre l'Exercice 8.4. D'autres stratégies
sont décrites par Huss-Lederman, Jacobson, Tsao, Turnbull et Johnson [Hus+96].

Bornes inférieures et bornes supérieures

Le problème des bornes inférieures est très difficile. En taille $n = 2$, il a été mon-
tré [HK71; Win71] que le résultat de Strassen est optimal du point de vue des multi-
plications : on ne peut pas effectuer le produit de deux matrices carrées de taille 2
en seulement 6 multiplications scalaires. Récemment Landsberg [Lan06] a obtenu
une généralisation de ce résultat, en utilisant des outils géométriques. Probert [Pro76]
a montré que la version de l'algorithme de Strassen due à Winograd est optimale
également vis-à-vis du nombre total d'opérations arithmétiques : il n'existe aucun
algorithme pour multiplier deux matrices $2 \times 2$ en 7 multiplications et seulement 14
additions et soustractions.

La taille $n = 2$ est la seule pour laquelle on connaît le nombre minimal de multi-
plications scalaires requises pour effectuer le produit de matrices carrées $n \times n$. En
taille $n = 3$, la meilleure borne inférieure connue à ce jour est 19 [Blä03], alors que
le meilleur algorithme commutatif utilise 22 multiplications [Mak86] et le meilleur
algorithme non commutatif utilise 23 multiplications [Lad76]. Pour une taille arbi-
traire, toutes les bornes inférieures connues sont seulement quadratiques en $n$ : les
meilleures, $\frac{5}{2}n^2 - 3n$ et $3n^2 - o(n^2)$, sont dues à Bläser [Blä01] et Landsberg [Lan14].
Si le corps des scalaires est $\mathbb{R}$ ou $\mathbb{C}$, et dans un modèle restreint de complexité où les
seules multiplications scalaires permises se font par des scalaires de module borné,
Raz [Raz03] a prouvé une borne inférieure en $cn^2 \log n$.

Matrices de petite taille

Il existe un grand nombre d'autres schémas à la *Strassen* pour multiplier des
matrices de petite taille. Un recensement des meilleurs algorithmes non commutatifs
connus pour le produit des matrices de petite taille a été effectué par Smith [Smi02].

Par exemple, Sýkora [Sýk77] a montré qu'on peut multiplier deux matrices $n \times n$
en $n^3 - (n-1)^2$ produits non commutatifs; cet algorithme généralise à la fois ceux de
Strassen et de Laderman. Des techniques plus sophistiquées ont permis à Pan [Pan78;
Pan80] de montrer que le produit de matrices $n \times n$ peut s'effectuer en $n^3/3+9n^2/2-n/3$
multiplications non commutatives. Chacun de ces schémas, appliqué récursivement,
donne un algorithme de produit matriciel sous-cubique. Par exemple, Pan [Pan80]
montre que l'on peut multiplier des matrices $48 \times 48$ en $47216$ multiplications non
commutatives, ce qui mène à $\omega_{mul} < 2{,}781$.
À la recherche de $\omega_{\text{mul}} = 2$

L'idée initiale de Strassen de combiner une bonne multiplication en taille fixée et la puissance de la récursivité pour améliorer l'exposant $\omega_{\text{mul}}$ n'a pas permis d'aller en dessous de 2,7. Les améliorations ultérieures reposent toutes sur une accumulation d'idées de plus en plus sophistiquées. Une première idée fut d'utiliser le produit des matrices rectangulaires. Si l'on dispose d'un algorithme pour multiplier deux matrices $(m \times n) \times (n \times p)$ en $\ell$ multiplications, alors on sait effectuer en $\ell$ multiplications chacun des produits en tailles $(n \times p) \times (p \times m)$ et $(p \times m) \times (m \times n)$ [HM73; Pan72]. Cet résultat est lié au théorème de transposition de Tellegen qui sera vu au Chapitre 12. Une décomposition par blocs permet alors d'obtenir un algorithme pour le produit de matrices carrées de taille $mnp$ en $\ell^3$ opérations, et d'obtenir ainsi la majoration $\omega_{\text{mul}} \le 3 \log(\ell)/\log(mnp)$. Ce résultat n'a permis d'améliorer la majoration de $\omega_{\text{mul}}$ que combiné avec d'autres concepts, comme les algorithmes approchés et le rang tensoriel limite (border rank en anglais) introduits par Bini et ses collaborateurs [Bin+79; Bin80]. L'article de Bini, Capovani, Romani et Lotti [Bin+79] donne un algorithme approché qui utilise 10 multiplications pour le produit matriciel en taille $(2 \times 2) \times (2 \times 3)$; cela implique $\omega_{\text{mul}} < \log_{12}(1000) < 2{,}77989$, mais les constantes dans les $O(\cdot)$ commencent à devenir astronomiques. D'autres idées pour obtenir des bornes de plus en plus fines sur $\omega_{\text{mul}}$ sont : l'utilisation des matrices creuses [Sch81] qui aboutit à l'inégalité asymptotique pour les sommes directes de matrices à trous et à la majoration $\omega_{\text{mul}} < 2{,}548$; la théorie de la déformation des modules sur des groupes algébriques linéaires conduisant à la méthode laser de Strassen [Str87] et à la borne $\omega_{\text{mul}} < 2{,}4785$, et qui est à la base de la borne $\omega_{\text{mul}} < 2{,}3729269$ due à Vassilevska Williams [Vas12; Vas14], et améliorant le précédent record dû à Coppersmith et Winograd [CW90]. La meilleure borne actuellement connue, $\omega_{\text{mul}} < 2{,}3728639$, est due à Le Gall [Le 14].

Théorie et pratique

Théoriquement, le principal problème ouvert est $\omega_{\text{mul}} = 2$. D'un point de vue pratique, il n'est pas forcément intéressant de comparer les algorithmes en fonction de leur exposant, mais plutôt de répondre à des questions du type : quelle est la formule la plus simple qui mène à un exposant réalisable $\theta < \log_2 7$ ? Quel est l'algorithme le plus rapide pour une taille donnée ?

L'idée des algorithmes mentionnés en Section 8.2 est due à Pan [Pan72; Pan81]. L'article de Laderman, Pan et Sha [LPS92] est le premier à avoir mis en avant leur aspect pratique; les derniers progrès basés sur cette idée sont dûs à Kaporin [Kap04; Kap99].

La recherche théorique pour améliorer les bornes sur $\omega_{\text{mul}}$ a demandé beaucoup d'ingéniosité (et de chance!), mais a mené à des algorithmes de plus en plus éloignés du calcul réel. Smith [Smi02] adopte un point de vue pragmatique et propose d'automatiser la recherche d'algorithmes rapides en petite taille. Il ramène la recherche d'un algorithme qui calcule le produit de deux matrices de types $(a \times b)$ et $(b \times c)$ en $m$ multiplications non commutatives à un problème d'optimisation dans un espace à $(ab + bc + ca)m$ dimensions. Cela lui permet de retrouver tous les algorithmes connus pour les petites tailles et d'en découvrir des nouveaux. Ses programmes retrouvent les formules de Strassen en quelques secondes et celles de Laderman en quelques minutes.
Un point de vue différent est celui développé par Landsberg [Lan08], qui montre que les questions ouvertes concernant la complexité de la multiplication matricielle se formulent convenablement en termes géométriques et théorie de la représentation.

Enfin, une autre piste actuellement explorée pour prouver (ou infirmer) que $\omega_{\text{mul}} = 2$ est celle ouverte par Cohn et Umans [Coh+05; CU03]. Ils développent une approche unifiée permettant de ramener la multiplication matricielle à l'étude de la transformée de Fourier discrète dans des algèbres de groupe associées à certains groupes finis, possédant de bonnes propriétés vis-à-vis de leurs représentations irréductibles. Ils retrouvent par ce biais $\omega_{\text{mul}} < 2,41$ et ramènent $\omega_{\text{mul}} = 2$ à des conjectures combinatoires et en théorie des groupes.

## Déterminant et permanent

Le permanent d'une matrice $A \in \mathcal{M}_n(\mathbb{K})$ est défini comme

$$ \operatorname{perm}(A) = \sum_{\sigma \in S_n} \prod_{i=1}^{n} a_{i, \sigma(i)}. $$

Malgré la ressemblance entre cette définition et celle de $\det(A)$, on ne connaît pas d'algorithmes de complexité polynomiale en $n$ calculant $\operatorname{perm}(A)$ sur un corps $\mathbb{K}$ de caractéristique différente de 2. Valiant [Val79] a prouvé que le calcul du permanent pour les matrices $n \times n$ à coefficients dans $\{0, 1\}$ est un problème NP-difficile.

## Bibliographie

AL04 ABDELJAOUED, J. et H. LOMBARDI (2004). *Méthodes matricielles : introduction à la complexité algébrique*. Vol. 42. Mathématiques & Applications. Springer-Verlag.

AL81 ALT, H. et J. van LEEUWEN (1981). « The complexity of basic complex operations ». In : *Computing. Archiv für Informatik und Numerik*, vol. 27, n°3, p. 205–215.

BCS97 BÜRGISSER, Peter, Michael CLAUSEN et M. Amin SHOKROLLAHI (1997). *Algebraic complexity theory*. Vol. 315. Grundlehren der Mathematischen Wissenschaften. Springer-Verlag.

Ber84 BERKOWITZ, Stuart J. (1984). « On computing the determinant in small parallel time using a small number of processors ». In : *Information Processing Letters*, vol. 18, p. 147–150.

BH74 BUNCH, James R. et John E. HOPCROFT (1974). « Triangular factorization and inversion by fast matrix multiplication ». In : *Mathematics of Computation*, vol. 28, p. 231–236.

Bin+79 BINI, D., M. CAPOVANI, F. ROMANI et G. LOTTI (1979). « $O(n^{2.7799})$ complexity for $n \times n$ approximate matrix multiplication ». In : *Information Processing Letters*, vol. 8, n°5, p. 234–235.

Bin80 BINI, D. (1980). « Relations between exact and approximate bilinear algorithms. Applications ». In : *Calcolo*, vol. 17, n°1, p. 87–97.

BK78 BRENT, R. P. et H. T. KUNG (1978). « Fast algorithms for manipulating formal power series ». In : *Journal of the ACM*, vol. 25, n°4, p. 581–595.
BKL91 BÜRGISSER, P., M. KARPINSKI et T. LICKTEIG (1991). « Some computational problems in linear algebra as hard as matrix multiplication ». In : *Computational Complexity*, vol. 1, n°2, p. 131-155.

Blä01 BLÄSER, Markus (2001). « A $\frac{5}{2}n^2$-lower bound for the multiplicative complexity of $n \times n$-matrix multiplication ». In : *STACS 2001 (Dresden)*. Vol. 2010. Lecture Notes in Computer Science. Springer-Verlag, p. 99-109.

Blä03 — (2003). « On the complexity of the multiplication of matrices of small formats ». In : *Journal of Complexity*, vol. 19, n°1, p. 43-60.

BM71 BORODIN, A. et I. MUNRO (1971). « Evaluating polynomials at many points ». In : *Information Processing Letters*, vol. 1, n°2, p. 66-68.

BS83 BAUR, W. et V. STRASSEN (1983). « The complexity of partial derivatives ». In : *Theoretical Computer Science*, vol. 22, p. 317-330.

Coh+05 COHN, H., R. KLEINBERG, B. SZEGEDY et C. UMANS (2005). « Group-theoretic algorithms for matrix multiplication ». In : *FOCS'05 : IEEE Conference on Foundations of Computer Science*. IEEE Computer Society, p. 379-388.

CU03 COHN, Henry et Christopher UMANS (2003). « A group-theoretic approach to fast matrix multiplication ». In : *FOCS'03 : IEEE Conference on Foundations of Computer Science*. Washington, DC, USA : IEEE Computer Society, p. 438.

CW90 COPPERSMITH, Don et Shmuel WINOGRAD (1990). « Matrix multiplication via arithmetic progressions ». In : *Journal of Symbolic Computation*, vol. 9, n°3, p. 251-280.

Dan37 DANILEVSKII, A. M (1937). « The numerical solution of the secular equation ». In : *Matematicheski Sbornik*, vol. 44, n°2. (in Russian), p. 169-171.

Fid72c FIDUCCIA, Charles M. (1972). « On obtaining upper bounds on the complexity of matrix multiplication ». In : *Complexity of computer computations (Proc. Sympos., IBM Thomas J. Watson Res. Center, Yorktown Heights, N.Y., 1972)*. New York : Plenum, p. 31-40.

Fis74 FISCHER, Patrick C. (1974). « Further schemes for combining matrix algorithms ». In : *Automata, languages and programming*. Vol. 14. Lecture Notes in Computer Science. Springer-Verlag, p. 428-436.

FP74 FISCHER, P. C. et R. L. PROBERT (1974). « Efficient procedures for using matrix algorithms ». In : *Automata, languages and programming*. Vol. 14. Lecture Notes in Computer Science. Springer-Verlag, p. 413-427.

Gie95 GIESBRECHT, Mark (1995). « Nearly optimal algorithms for canonical matrix forms ». In : *SIAM Journal on Computing*, vol. 24, n°5, p. 948-969.

Har72 HARTER, Richard (1972). « The optimality of Winograd's formula ». In : *Communications of the ACM*, vol. 15, n°5, p. 352.

HK71 HOPCROFT, J. E. et L. R. KERR (1971). « On minimizing the number of multiplications necessary for matrix multiplication ». In : *SIAM Journal on Applied Mathematics*, vol. 20, p. 30-36.

HM73 HOPCROFT, J. et J. MUSINSKI (1973). « Duality applied to the complexity of matrix multiplication and other bilinear forms ». In : *SIAM Journal on Computing*, vol. 2, p. 159-173.

Hus+96 HUSS-LEDERMAN, S., E. M. JACOBSON, A. TSAO, T. TURNBULL et J. R. JOHNSON (1996). « Implementation of Strassen's algorithm for matrix multiplication ».
tion». In: Supercomputing'96. 32 pp. Pittsburgh, PA: IEEE Computer Society.

JaJ79 Ja'Ja', Joseph (1979). « On the complexity of bilinear forms with commutativity ». In: STOC'79: ACM Symposium on Theory of Computing. Atlanta, Georgia, United States: ACM, p. 197–208.

Kap04 KAPORIN, I. (2004). « The aggregation and cancellation techniques as a practical tool for faster matrix multiplication ». In: Theoretical Computer Science, vol. 315, n°2-3, p. 469–510.

Kap99 KAPORIN, Igor (1999). « A practical algorithm for faster matrix multiplication ». In: Numerical Linear Algebra with Applications, vol. 6, n°8, p. 687–700.

Kel85 KELLER-GEHRIG, Walter (1985). « Fast algorithms for the characteristic polynomial ». In: Theoretical Computer Science, vol. 36, n°2-3, p. 309–317.

KK65 KLYUYEV, V. V. et N. I. KOKOVKIN-SHCHERBAK (1965). « On the minimization of the number of arithmetic operations for the solution of linear algebraic systems of equations ». In: USSR Computational Mathematics and Mathematical Physics, vol. 5, p. 25–43.

Kry31 KRYLOV, A. N. (1931). « On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined ». In: Izvestiya Akademii Nauk SSSR, vol. 7, n°4. (in Russian), p. 491–539.

Lad76 LADERMAN, Julian D. (1976). « A noncommutative algorithm for multiplying 3×3 matrices using 23 multiplications ». In: Bulletin of the American Mathematical Society, vol. 82, n°1, p. 126–128.

Lan06 LANDSBERG, J. M. (2006). « The border rank of the multiplication of 2 × 2 matrices is seven ». In: Journal of the American Mathematical Society, vol. 19, n°2, p. 447–459.

Lan08 — (2008). « Geometry and the complexity of matrix multiplication ». In: Bulletin of the American Mathematical Society, vol. 45, n°2, p. 247–284.

Lan14 — (2014). « New lower bounds for the rank of matrix multiplication ». In: SIAM Journal on Computing, vol. 43, n°1, p. 144–149.

Le 14 LE GALL, François (2014). « Powers of tensors and fast matrix multiplication ». In: ISSAC'14: International Symposium on Symbolic and Algebraic Computation. Éd. par Katsusuke NABESHIMA. New York, NY, USA: ACM, p. 296–303.

Le 40 LE VERRIER, Urbain (1840). « Sur les variations séculaires des éléments elliptiques des sept planètes principales: Mercure, Vénus, la Terre, Mars, Jupiter, Saturne et Uranus ». In: Journal de Mathématiques Pures et Appliquées, vol. 1 (5), p. 220–254.

Lic87 LICKTEIG, Thomas (1987). « The computational complexity of division in quadratic extension fields ». In: SIAM Journal on Computing, vol. 16, n°2, p. 278–311.

LPS92 LADERMAN, Julian, Victor PAN et Xuan He SHA (1992). « On practical algorithms for accelerated matrix multiplication ». In: Linear Algebra and its Applications, vol. 162/164, p. 557–588.
Mak86 MAKAROV, O. M. (1986). « An algorithm for multiplying 3 × 3 matrices ». In : *USSR Computational Mathematics and Mathematical Physics*, vol. 26, n°1, p. 179–180.

Mun73 MUNRO, I. (1973). « Problems related to matrix multiplication ». In : *Proceedings Courant Institute Symposium on Computational Complexity*, October 1971. Éd. par R. RUSTIN. New York : Algorithmics Press, p. 137–151.

Pan72 PAN, V. (1972). « Computation schemes for a product of matrices and for the inverse matrix ». In : *Akademiya Nauk SSSR i Moskovskoe Matematicheskoe Obshchestvo. Uspekhi Matematicheskikh Nauk*, vol. 27, n°5(167), p. 249–250.

Pan78 PAN, V. Ya. (1978). « Strassen's algorithm is not optimal. Trilinear technique of aggregating, uniting and canceling for constructing fast algorithms for matrix operations ». In : *SFCS '78*. Washington, DC, USA : IEEE Computer Society, p. 166–176.

Pan80 — (1980). « New fast algorithms for matrix operations ». In : *SIAM Journal on Computing*, vol. 9, n°2, p. 321–342.

Pan81 — (1981). « New combinations of methods for the acceleration of matrix multiplication ». In : *Computers & Mathematics with Applications. An International Journal*, vol. 7, n°1, p. 73–125.

Pan84a PAN, Victor (1984). « How can we speed up matrix multiplication? » In : *SIAM Review*, vol. 26, n°3, p. 393–415.

Pan84b — (1984). *How to multiply matrices faster*. Vol. 179. Lecture Notes in Computer Science. Springer-Verlag.

Pat75 PATERSON, Michael S. (1975). « Complexity of monotone networks for Boolean matrix product ». In : *Theoretical Computer Science*, vol. 1, n°1, p. 13–20.

Pro76 PROBERT, Robert L. (1976). « On the additive complexity of matrix multiplication ». In : *SIAM Journal on Computing*, vol. 5, n°2, p. 187–203.

PS07 PERNET, Clément et Arne STORJOHANN (2007). « Faster algorithms for the characteristic polynomial ». In : *ISSAC'07 : International Symposium on Symbolic and Algebraic Computation*. Éd. par Dongming WANG. ACM Press, p. 307–314.

PS73 PATERSON, M. S. et L. J. STOCKMEYER (1973). « On the number of nonscalar multiplications necessary to evaluate polynomials ». In : *SIAM Journal on Computing*, vol. 2, n°1, p. 60–66.

Raz03 RAZ, Ran (2003). « On the complexity of matrix product ». In : *SIAM Journal on Computing*, vol. 32, n°5, p. 1356–1369.

Sch73 SCHÖNHAGE, A. (1972/73). « Unitäre Transformationen grosser Matrizen ». In : *Numerische Mathematik*, vol. 20, p. 409–417.

Sch81 — (1981). « Partial and total matrix multiplication ». In : *SIAM Journal on Computing*, vol. 10, n°3, p. 434–455.

Smi02 SMITH, Warren D. (2002). « Fast matrix multiplication formulae : report of the prospectors ». Preprint. http://www.math.temple.edu/~wds/prospector.pdf.

Smi62 SMITH, R. L. (1962). « Algorithm 116 : Complex division ». In : *Communications of the ACM*, n°5, p. 435.
Sto01 STORJOHANN, Arne (2001). « Deterministic computation of the Frobenius form ». In : FOCS'01 : IEEE Symposium on Foundations of Computer Science. Extended abstract. Las Vegas, NV : IEEE Computer Society, p. 368–377.

Str69 STRASSEN, V. (1969). « Gaussian elimination is not optimal ». In : Nume-rische Mathematik, vol. 13, p. 354–356.

Str87 — (1987). « Relative bilinear complexity and matrix multiplica- tion ». In : Journal für die Reine und Angewandte Mathematik, vol. 375/376, p. 406–443.

Str90 STRASSEN, Volker (1990). « Algebraic complexity theory ». In : Handbook of theoretical computer science, Vol. A. Amsterdam : Elsevier, p. 633–672.

Sýk77 SÝKORA, Ondrej (1977). « A fast non-commutative algorithm for matrix multiplication ». In : Mathematical foundations of computer science (Proc. Sixth Sympos., Tatranská Lomnica, 1977). Vol. 53. Lecture Notes in Computer Science. Springer-Verlag, p. 504–512.

Val79 VALIANT, L. G. (1979). « The complexity of computing the permanent ». In : Theoretical Computer Science, vol. 8, n°2, p. 189–201.

Vas12 VASSILEVSKA WILLIAMS, Virginia (2012). « Multiplying matrices faster than Coppersmith–Winograd ». In : STOC'12 : ACM Symposium on Theory of Computing. New York, NY : ACM, p. 887–898.

Vas14 — (2014). « Multiplying matrices in O(n²³⁷³) time ». http://theory.stanford.edu/~virgi/matrixmult-f.pdf.

Wak70 WAKSMAN, Abraham (1970). « On Winograd's algorithm for inner products ». In : IEEE Transactions on Computers, vol. C-19, n°4, p. 360–361.

Win67 WINOGRAD, S. (1967). « On the number of multiplications required to compute certain functions ». In : Proceedings of the National Academy of Sciences of the United States of America, vol. 58, p. 1840–1842.

Win68 — (1968). « A new algorithm for inner-product ». In : IEEE Transactions on Computers, vol. 17, p. 693–694.

Win71 — (1971). « On multiplication of 2 × 2 matrices ». In : Linear Algebra and Applications, vol. 4, p. 381–388.
# 9. Algèbre linéaire creuse : algorithme de Wiedemann

## Résumé

Les matrices creuses sont les matrices contenant beaucoup d'éléments nuls. L'algorithme de Wiedemann est une méthode itérative pour résoudre des systèmes linéaires représentés par des matrices creuses. Il ramène la résolution au calcul du polynôme minimal, lui-même reposant sur la reconnaissance d'une suite récurrente à coefficients constants.

## 9.1 Introduction

Dans ce chapitre, on aborde des questions radicalement différentes de celles du Chapitre 8 : on ne cherche plus à effectuer les produits, inversions, ... de matrices quelconques, mais à manipuler des matrices creuses.

On ne commencera pas par définir précisément ce qu'est une matrice creuse. L'aproche consiste à donner des algorithmes dont la complexité s'exprime en fonction du nombre d'éléments non nuls des matrices en entrée. De la complexité des algorithmes « creux » va découler une borne, typiquement de l'ordre $O(n)$ pour des matrices de taille $n \times n$, sur le nombre de coefficients non nuls pour que le changement de modèle soit pertinent.

Dans tout ce qui suit, on considère donc une matrice $A$ de taille $n \times n$ à coefficients dans un corps $\mathbb{K}$ et vérifiant les hypothèses suivantes :

- A est inversible,

- A contient $s$ éléments non nuls.

Le produit de $A$ par un vecteur peut donc s'effectuer en $O(s)$ opérations dans $\mathbb{K}$.

On va décrire un algorithme dû à Wiedemann qui permet de calculer l'unique solution du système $Ax = y$ avec une complexité en $O(ns)$ ; l'algorithme original
est plus général que celui présenté ici, en ce qu'il permet de traiter les matrices rectangulaires ou carrées et non inversibles, mais il se ramène au cas carré inversible.

Remarquons que si $s$ est de l'ordre de $n^2$, caractérisant donc des matrices plutôt denses, on retombe dans une complexité de l'ordre de $O(n^3)$. Le cas intéressant est celui où $s$ est de l'ordre de $n$, auquel cas l'algorithme est quadratique en $n$ : pour simplifier, on retiendra que l'exposant de l'algèbre linéaire creuse est 2 dans les bons cas.

## 9.2 Polynôme minimal et résolution de systèmes

L'algorithme de Wiedemann passe par le calcul du polynôme minimal de A. Rap-
pelons sa définition. Il est possible d'associer à tout polynôme en une variable
$P = \sum_i p_i X^i$ de $\mathbb{K}[X]$ la matrice $P(A) = \sum_i p_i A^i$. Le théorème de Cayley-Hamilton
affirme que le polynôme caractéristique $\chi_A$ de A annule la matrice A, c'est-à-dire que
$\chi_A(A) = 0$.

**Définition 9.1** Le polynôme minimal de $A$ est le polynôme unitaire de plus petit degré annulant $A$.

Soit $\mu_A(X)$ le polynôme minimal de $A$, supposé connu. Puisqu'on a supposé $A$ inversible, le terme constant de $\mu_A(X)$ n'est pas nul : il divise le terme constant du polynôme caractéristique, qui n'est lui-même pas nul.

L'égalité $\mu_A(A) = 0$ se récrit sous la forme

$$A^{-1} = \frac{-1}{p_0} (A^{m-1} + p_{m-1} A^{m-2} + \cdots + p_1 I_n).$$

Pour résoudre le système $Ax = b$, $(b \in \mathbb{K}^n)$, on va calculer $x = A^{-1}b$, c'est-à-dire

$$\frac{-1}{p_0} (A^{m-1} b + p_{m-1} A^{m-2} b + \cdots + p_1 b).$$

Ainsi, il suffit de calculer les itérés $b, Ab, \dots, A^{m-1}b$, puis d'additionner les vecteurs $p_1b, p_2Ab, \dots, A^{m-1}b$, pour retrouver $x$. L'analyse de complexité est immédiate :

— Chacun des $A^i b$ s'obtient à partir de $A^{i-1} b$ en $O(s)$ opérations.

— Chacun des produits par $p_i$, et chaque addition, se fait en $O(n)$ opérations.

Remarquer que $n \le s$ (hypothèse d'inversibilité de $A$).

Au total, on a donc $O(ns)$ opérations à faire pour résoudre le système, si on connaît le polynôme minimal de $A$.

## 9.3 Calcul du polynôme minimal

Écrivons le polynôme minimal de $A$ sous la forme

$$\mu_A(X) = X^m + p_{m-1} X^{m-1} + \cdots + p_0.$$

Alors la suite des puissances de $A$ vérifie la récurrence linéaire

$$A^{k+m} + p_{m-1} A^{k+m-1} + \cdots + p_0 A^k = 0,$$
**Entrée** Une matrice creuse $A \in \mathcal{M}_n(\mathbb{K})$.

**Sortie** Son polynôme minimal $\mu_A(X)$.

1. Choisir des vecteurs $u$ et $v$ aléatoires dans $\mathbb{K}^n$.

2. Pour $0 \le i \le 2n$, calculer les vecteurs $v_i = A^i v$, puis les scalaires $N_i = {}^t u v_i$.

3. Renvoyer le polynôme minimal de la suite $(N_i)$.

**Algorithme 9.1** – Algorithme de Wiedemann pour calculer le polynôme minimal d'une matrice creuse.

et ne vérifie aucune récurrence d'ordre plus petit.

Pour tous vecteurs $u$ et $v$, la suite (de nombres) $N_i := {}^t u A^i v$ vérifie donc

$$N_{k+m} + p_{m-1}N_{k+m-1} + \dots + p_0N_k = 0.$$

Si $u$ et $v$ sont mal choisis (nuls, par exemple), la suite $N_i$ vérifie une récurrence d'ordre plus petit que $m$. La proposition suivante montre qu'en choisissant $u$ et $v$ au hasard, on tombe vraisemblablement sur une suite $N_i$ qui ne satisfait pas de récurrence d'ordre plus petit.

**Proposition 9.1** Il existe un polynôme $D \in \mathbb{K}[U_1, \dots, U_n, V_1, \dots, V_n]$, non nul et de degré au plus $2n$, tel que si $D(u, v)$ est non nul, la suite $N_i$ associée à $u$ et $v$ ne satisfait pas de récurrence d'ordre plus petit que $m$.

L'idée est de choisir $u$ et $v$ au hasard. Si on n'a pas de chance, $D(u, v)$ est nul; sinon, la suite des $N_i$ associée à $u$ et $v$ permet de retrouver $\mu_A(X)$ grâce à l'algorithme de Berlekamp-Massey décrit en page 143 (attention, on ne connaît pas explicitement le polynôme $D$, mais son existence assure que la plupart des choix de $u$ et $v$ sont « chanceux »).

L'Algorithme 9.1 résume cette approche probabiliste. L'analyse de probabilité est facile à faire en utilisant les résultats de l'Exercice 9.1 ci-dessous.

La complexité de l'étape (2) est de $O(n)$ produits matrice-vecteur ou vecteur-vecteur, chacun prenant au pire $O(s)$ opérations; au total on obtient donc $O(ns)$ opérations pour cette phase. L'algorithme d'approximants de Padé est négligeable, puisqu'il ne demande que $O(n^2) = O(sn)$ opérations.

## 9.4 Calcul du déterminant

L'Algorithme 9.1 permet de détecter en cours de calcul l'inversibilité de la matrice $A$. En effet, si l'approximant de Padé calculé par l'algorithme de Berlekamp-Massey (page 143) a un dénominateur divisible par $X$, cela veut dire que $\det(A) = 0$.

Il est possible de calculer en bonne complexité le déterminant de la matrice creuse $A$, en exploitant le résultat suivant d'algèbre linéaire.
**Lemme 9.2** Si tous les mineurs principaux de $A \in M_n(\mathbb{K})$ sont non nuls, alors les polynômes minimal et caractéristique de la matrice $B = A \cdot \text{Diag}(X_1, \dots, X_n)$ coïncident.

L'algorithme qui s'en déduit est le suivant :

— Choisir une matrice de permutations aléatoire $P$, pour assurer que tous les mineurs principaux de la matrice $AP$ sont inversibles;

— choisir des éléments $x_1, \dots, x_n$ aléatoirement dans $\mathbb{K}$;

— calculer le polynôme minimal $\mu_B$ de $B = A \cdot P \cdot \text{Diag}(x_1, \dots, x_n)$;

— renvoyer $\mu_B(0)/(x_1 \cdots x_n) \cdot \text{sgn}(P)$.

À nouveau, la complexité de cet algorithme est quadratique en la taille $n$ de $A$.

## 9.5 Calcul du rang

Le rang maximal de la matrice creuse $A$ est détecté par les algorithmes précédents. Quitte à multiplier $A$ par une matrice de permutations aléatoires, on peut donc supposer que son rang $r$ est strictement inférieur à $n$, et que ses mineurs principaux $A_i$, pour $1 \le i \le r$, sont tous non nuls.

Sous ces hypothèses, il est possible de montrer que, pour un choix aléatoire d'éléments $x_1, \dots, x_n$ dans $\mathbb{K}$, le rang $r$ est égal à $\text{deg}(\mu_{AD}) - 1$, où $D$ est la matrice diagonale $D = \text{Diag}(x_1, \dots, x_n)$. L'algorithme qui s'ensuit calcule le polynôme minimal de la matrice creuse AD, et en déduit le rang de $A$, le tout pour une complexité quadratique en $n$.

## Exercices

**Exercice 9.1** Soit $A$ une matrice dans $M_n(\mathbb{K})$, soit $b$ un vecteur de $\mathbb{K}^n \setminus \{0\}$ et soit $f$ le polynôme minimal de la suite $(A^i \cdot b)_{i \ge 0}$.

Le but de l'exercice est d'estimer la probabilité $P$ que $f$ coïncide avec le polynôme minimal de la suite $(^t u \cdot A^i \cdot b)_{i \ge 0}$ lorsque $u$ est un vecteur de $\mathbb{K}^n$ dont les coordonnées sont choisies aléatoirement au hasard dans un sous-ensemble fini $U$ de $\mathbb{K}$.

1. Montrer qu'il existe une application $\psi : \mathbb{K}^n \to \mathbb{A} = \mathbb{K}[X]/(f)$, $\mathbb{K}$-linéaire et surjective, telle que pour tout $u \in \mathbb{K}^n$ on ait

$$f \text{ est le polynôme minimal de } (^t u \cdot A^i \cdot b)_{i \ge 0} \iff \psi(u) \text{ est inversible dans } \mathbb{A}.$$

[Indication : l'application $\phi : \mathbb{K}^n \to \mathbb{K}^N$ définie par $\phi(u) = (^t u \cdot A^i \cdot b)_{i \ge 0}$ induit une application linéaire surjective $\mathbb{K}^n \to M_f$, où $M_f$ est l'espace vectoriel des suites de $\mathbb{K}^N$ admettant $f$ comme polynôme annulateur. Par ailleurs, $\mathbb{A}$ et $M_f$ sont isomorphes en tant qu'espaces vectoriels.]

2. Soit $d$ le degré de $f$. Montrer qu'il existe un polynôme non identiquement nul $R \in \mathbb{K}[X_1, \dots, X_n]$ de degré total au plus $d$ tel que pour tout vecteur $u = (u_1, \dots, u_n) \in \mathbb{K}^n$ on ait

$$\psi(u) \text{ est inversible dans } \mathbb{A} \text{ si et seulement si } R(u_1, \dots, u_n) \neq 0.$$
[Indication : utiliser un résultant.]

3. Montrer que R admet au plus $d \cdot \text{card}(U)^{n-1}$ racines dans $U^n$. En déduire que la probabilité qu'un élément de $U^n$ dont les coordonnées sont choisies aléatoirement au hasard dans U soit racine de R est bornée par $d/\text{card}(U)$.

4. Conclure que la probabilité $\mathcal{P}$ vérifie

$$ \mathcal{P} \geq 1 - \frac{d}{\text{card}(U)}. $$

■

## Notes

L'algorithme de Wiedemann est décrit dans son article de 1986 [Wie86], qui contient également l'algorithme esquissé en Section 9.4. L'algorithme en Section 9.5 est dû à Kaltofen et Saunders [KS91]. Des versions « par blocs » de l'algorithme de Wiedemann ont été proposées par Coppersmith [Cop94], Villard [Vil97], Thomé [Tho02], et Turner [Tur06]. D'autres généralisations se trouvent dans les travaux de Giesbrecht, Lobo et Saunders [Gie01; GLS98], et de Mulders [Mul04].

Les questions (2) et (3) de l'Exercice 9.1 forment le cœur du « lemme de Schwartz-Zippel » [DL78; Sch80b; Zip79]. Nous en verrons une preuve en page 410 du Chapitre 22 sur la factorisation des polynômes à plusieurs variables.

## Bibliographie

Cop94 COPPERSMITH, Don (1994). « Solving homogeneous linear equations over GF(2) via block Wiedemann algorithm ». In : Mathematics of Computation, vol. 62, n°205, p. 333-350.

DL78 DEMILLO, Richard A. et Richard J. LIPTON (1978). « A probabilistic remark on algebraic program testing ». In : Information Processing Letters, vol. 7, n°4, p. 193-195.

Gie01 GIESBRECHT, Mark (2001). « Fast computation of the Smith form of a sparse integer matrix ». In : Computational Complexity, vol. 10, n°1, p. 41-69.

GLS98 GIESBRECHT, M., A. LOBO et B. D. SAUNDERS (1998). « Certifying inconsistency of sparse linear systems ». In : ISSAC'98 : International Symposium on Symbolic and Algebraic Computation. New York : ACM Press, p. 113-119.

KS91 KALTOFEN, Erich et B. David SAUNDERS (1991). « On Wiedemann's method of solving sparse linear systems ». In : Applied algebra, algebraic algorithms and error-correcting codes (New Orleans, LA, 1991). Vol. 539. Lecture Notes in Computer Science. Springer-Verlag, p. 29-38.

Mul04 MULDERS, Thom (2004). « Certified sparse linear system solving ». In : Journal of Symbolic Computation, vol. 38, n°5, p. 1343-1373.

Sch80b SCHWARTZ, J. T. (1980). « Fast probabilistic algorithms for verification of polynomial identities ». In : Journal of the Association for Computing Machinery, vol. 27, n°4, p. 701-717.
Tho02 THOMÉ, E. (2002). « Subquadratic computation of vector generating polynomials and improvement of the block Wiedemann algorithm ». In : *Journal of Symbolic Computation*, vol. 33, n°5, p. 757-775.

Tur06 TURNER, William J. (2006). « A block Wiedemann rank algorithm ». In : *ISSAC'06 : International Symposium on Symbolic and Algebraic Computation*. New York : ACM Press, p. 332-339.

Vil97 VILLARD, G. (1997). « Further analysis of Coppersmith's block Wiedemann algorithm for the solution of sparse linear systems ». In : *ISSAC'97 : International Symposium on Symbolic and Algebraic Computation*. Kihei, Maui, Hawaii, United States : ACM Press, p. 32-39.

Wie86 WIEDEMANN, D. (1986). « Solving sparse linear equations over finite fields ». In : *IEEE Transactions on Information Theory*, vol. IT-32, p. 54-62.

Zip79 ZIPPEL, Richard (1979). « Probabilistic algorithms for sparse polynomials ». In : *Symbolic and Algebraic Computation. Eurosam '79, An International Symposium on Symbolic and Algebraic Manipulation, Marseille, France, June 1979*. Vol. 72. Lecture Notes in Computer Science. Springer-Verlag, p. 216-226.
# 10. Algèbre linéaire structurée

## Résumé

Les matrices possédant une structure spéciale sont omniprésentes en calcul formel. Ce sont des matrices dont les éléments jouissent d'une certaine répétition, ou satisfont à certaines relations. Plus formellement, une matrice $n \times n$ structurée est typiquement définie par $O(n)$ éléments, au lieu de $n^2$ pour une matrice dense, et peut être multipliée par un vecteur en $\mathcal{O}(n)$ opérations arithmétiques, au lieu de $O(n^2)$ opérations pour une matrice dense. Ce chapitre présente une algorithmique unifiée pour manipuler des matrices structurées denses, comme les matrices de Toeplitz, de Hankel, de Vandermonde et de Sylvester. Les calculs avec les matrices de ces classes sont liés aux calculs avec les polynômes, ce qui permet l'emploi des techniques de multiplication polynomiale rapide pour accélérer leur manipulation. Par exemple, on peut résoudre un système linéaire défini par une matrice structurée inversible en $\mathcal{O}(n)$ opérations.

## 10.1 Introduction

Nous avons vu au Chapitre 8 qu'il est possible de manipuler les matrices denses de taille $n \times n$ à coefficients dans un corps $\mathbb{K}$ en $MM(n) = O(n^\theta)$ opérations dans $\mathbb{K}$, où $\theta$ est un réel compris entre 2 et 3. Ici, par manipuler, on entend multiplier, inverser, calculer le déterminant, ou encore résoudre un système linéaire.

Dans certaines situations, les matrices qu'on est amené à manipuler présentent une structure que ces algorithmes généraux ne savent pas détecter et exploiter. Voici quelques exemples de matrices qui possèdent une structure.
**Définition 10.1** Une matrice $A \in \mathcal{M}_n(\mathbb{K})$ est dite de type Toeplitz (ou simplement de Toeplitz) si elle est invariante le long des diagonales, c'est-à-dire si ses éléments $a_{i,j}$ vérifient $a_{i,j} = a_{i+k,j+k}$ pour tout $k$.

Une telle matrice est complètement définie par sa première ligne et sa première colonne.

**Définition 10.2** Une matrice $A \in \mathcal{M}_n(\mathbb{K})$ est dite de Hankel si elle est invariante le long des anti-diagonales, c'est-à-dire si ses éléments $a_{i,j}$ vérifient $a_{i,j} = a_{i-k,j+k}$ pour tout $k$.

Une telle matrice est complètement définie par sa première ligne et sa dernière colonne.

**Exemple 10.1** La matrice de la multiplication polynomiale en degré fixé (dans les bases canoniques) par un polynôme fixé est de Toeplitz. Par exemple, la multiplication d'un polynôme de degré au plus 2 par le polynôme fixé $a_0 + a_1 X + a_2 X^2$ de $\mathbb{K}[X]$ se traduit matriciellement par l'égalité

$$ \begin{pmatrix} a_0 & 0 & 0 \\ a_1 & a_0 & 0 \\ a_2 & a_1 & a_0 \\ 0 & a_2 & a_1 \\ 0 & 0 & a_2 \end{pmatrix} \times \begin{pmatrix} b_0 \\ b_1 \\ b_2 \end{pmatrix} = \begin{pmatrix} a_0 b_0 \\ a_0 b_1 + a_1 b_0 \\ a_0 b_2 + a_1 b_1 + a_2 b_0 \\ a_1 b_2 + a_2 b_1 \\ a_2 b_2 \end{pmatrix} $$

Dans cet exemple, la matrice de Toeplitz est d'une forme particulière, appelée bande. En fait, toute matrice de Toeplitz peut être vue comme sous-matrice d'une matrice de Toeplitz bande, et cela entraîne le résultat suivant.

**Lemme 10.1** Le produit d'une matrice de Toeplitz (ou de Hankel) de $\mathcal{M}_n(\mathbb{K})$ par un vecteur de $\mathbb{K}^n$ peut s'effectuer en $O(M(n))$ opérations.

*Démonstration.* Pour tout $0 \le i \le n-1$, l'élément $c_i$ du produit matrice-vecteur

$$ \begin{pmatrix} a_{n-1} & \cdots & a_0 \\ \vdots & \ddots & \vdots \\ a_{2n-2} & \cdots & a_{n-1} \end{pmatrix} \times \begin{pmatrix} b_0 \\ \vdots \\ b_{n-1} \end{pmatrix} = \begin{pmatrix} c_0 \\ \vdots \\ c_{n-1} \end{pmatrix} $$

est le coefficient de $X^{n-1+i}$ dans le produit de polynômes

$$ (a_0 + \cdots + a_{2n-2}X^{2n-2}) \times (b_0 + \cdots + b_{n-1}X^{n-1}). \quad (10.1) $$

La preuve est tout à fait similaire pour une matrice de Hankel.

<figure><img src="image_28.png" /></figure>
R

La preuve précédente montre que $2M(n) + O(n)$ opérations suffisent pour multiplier une matrice de Toeplitz ou de Hankel de taille $n$ par un vecteur. Cette borne peut être réduite à $M(n) + O(n)$, en observant que l'extraction de la partie médiane du produit (10.1) est suffisante, et en employant des algorithmes pour le produit médian, évoqués au Chapitre 12.

**Exemple 10.2** Une matrice de Sylvester est la concaténation de deux matrices Toeplitz bande, cf. 6.2.

**Définition 10.3** Une matrice $A = (a_{i,j})_{i,j=0}^{n-1}$ de $M_n(\mathbb{K})$ est dite de Vandermonde si ses éléments s'écrivent $a_{i,j} = a_i^j$ pour $a_0, \dots, a_{n-1} \in \mathbb{K}$.

**Définition 10.4** Une matrice $A = (a_{i,j})_{i,j=0}^{n-1}$ de $M_n(\mathbb{K})$ est dite de Cauchy si ses éléments s'écrivent $a_{i,j} = 1/(a_i - b_j)$ pour $a_i, b_j \in \mathbb{K}$ avec $a_i \neq b_j$ pour tous $i, j$.

Il y a un certain nombre de points en commun entre ces exemples : la matrice est représentable par $O(n)$ éléments; le produit matrice-vecteur peut s'effectuer plus rapidement que dans le cas générique, en $O(M(n) \log n)$ opérations au lieu de $O(n^2)$; le produit par une matrice quelconque peut s'effectuer en complexité $O(nM(n) \log n)$ — quasi-optimale en la taille de la sortie pour une multiplication polynomiale à base de FFT. En effet, dans chaque cas, le produit matrice-vecteur $Av$ admet une interprétation analytique :

- multiplication polynomiale dans les cas Toeplitz, Hankel et Sylvester;

- évaluation polynomiale multipoint dans le cas d'une matrice de Vandermonde (Chapitre 5);

- évaluation multipoint de fractions rationnelles de la forme $\sum_{j=0}^{n-1} c_j/(X - b_j)$ dans le cas d'une matrice de Cauchy.

**Exercice 10.1** Montrer que le produit d'une matrice de Cauchy par un vecteur peut s'effectuer en $O(M(n) \log n)$ opérations dans $\mathbb{K}$. ■

Au vu de ces exemples, on pourrait donc être tenté de définir comme *structurée* une matrice telle que son produit par un vecteur peut s'effectuer en $\tilde{O}(n)$ opérations. La question qui se pose naturellement est : peut-on exploiter cette définition de la structure de $A$ pour résoudre le système $Ax = b$ ? Un premier constat positif est que, dans chacun des exemples précédents, la résolution admet aussi une interprétation analytique :

i. devinette de récurrences à coefficients constants, si $A$ est de Toeplitz ou de Hankel;

ii. interpolation polynomiale, si $A$ est de Vandermonde;

iii. interpolation de fractions rationnelles, si $A$ est de Cauchy.

En effet, si la suite $(a_n)$ d'éléments de $\mathbb{K}$ vérifie une récurrence (inconnue) à coefficients constants de la forme

$$a_{n+d} = p_{d-1}a_{n+d-1} + \cdots + p_0a_n, \quad n \ge 0,$$
alors trouver les coefficients $p_i$ de cette récurrence revient à résoudre le système de Hankel

$$ \begin{pmatrix} a_0 & a_1 & \cdots & a_{d-1} \\ a_1 & a_2 & \cdots & a_d \\ \vdots & \vdots & \ddots & \vdots \\ a_{d-1} & a_d & \cdots & a_{2d-2} \end{pmatrix} \times \begin{pmatrix} p_0 \\ p_1 \\ \vdots \\ p_{d-1} \end{pmatrix} = \begin{pmatrix} a_d \\ a_{d+1} \\ \vdots \\ a_{2d-1} \end{pmatrix}. $$

Un autre point encourageant est que, pour ces trois opérations, on dispose d'algorithmes de complexité quasi-optimale $O(M(n) \log n)$, présentés au Chapitre 5 pour (ii), et au Chapitre 7 pour (i) et (iii).

Il y a cependant quelques points négatifs : un premier est que si A est une matrice inversible telle que l'application linéaire $v \mapsto Av$ se calcule en L opérations, cela n'implique pas l'existence d'un algorithme de complexité $O(L)$ pour l'application $v \mapsto A^{-1}v$. En d'autres termes, on manque d'un principe d'inversion analogue au principe de transposition présenté au Chapitre 12. Par exemple, pour les matrices creuses, le meilleur algorithme de résolution, dû à Wiedemann, est de complexité quadratique en $n$ (Chapitre 9). Un second point négatif est que les classes vues plus haut ne sont pas stables par inversion : l'inverse d'une matrice de Toeplitz (resp. de Vandermonde) n'est pas de Toeplitz (resp. de Vandermonde).

On a donc besoin de rajouter des hypothèses dans la définition de la *bonne notion* de *matrice structurée*. Cette définition exploite la généralisation du caractère invariant par diagonale d'une matrice de Toeplitz, et vient de l'observation simple suivante : si A est de Toeplitz, alors la matrice

$$ \phi(A) = A - (A \text{ décalée de } 1 \text{ vers le bas et de } 1 \text{ vers la droite}) $$

admet une sous-matrice carrée nulle d'ordre un de moins; elle est donc de rang au plus 2 (on convient qu'on remplit avec des zéros les premières ligne et colonne de la matrice décalée). On dit alors que $\phi$ est un *opérateur de déplacement* et que le *rang de déplacement par rapport à* $\phi$ de A est au plus 2. On peut donc représenter $\phi(A)$ sous forme compacte, comme un produit $G \cdot {}^t H$, avec G et H des matrices rectangulaires de taille $n \times 2$.

**Exemple 10.3** Si A est la matrice $3 \times 3$

$$ A = \begin{pmatrix} c & d & e \\ b & c & d \\ a & b & c \end{pmatrix}, $$

avec $d \neq 0$, alors $\phi(A)$ s'écrit

$$ \phi(A) = \begin{pmatrix} c & d & e \\ b & 0 & 0 \\ a & 0 & 0 \end{pmatrix} = \begin{pmatrix} c & d \\ b & 0 \\ a & 0 \end{pmatrix} \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & e/d \end{pmatrix}. $$
**Définition 10.5** Les matrices (G, H) de $M_{n,2}(\mathbb{K})$ telles que $\phi(A) = G \cdot {}^t H$ sont appelées générateurs de déplacement pour la matrice de Toeplitz A.

Ces définitions s'étendent à une matrice quelconque et conduisent au concept de matrice *quasi-Toeplitz*.

**Définition 10.6** On appelle plus généralement *opérateur de déplacement* $\phi_+$ l'application $A \mapsto A - Z \cdot A \cdot {}^t Z$, où A est une matrice quelconque et Z la matrice définie par

$$Z = \begin{pmatrix} 0 & 0 & \dots & 0 \\ 1 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 1 & 0 \end{pmatrix}$$

ZA est la matrice A décalée de 1 ligne vers le bas, et $A \cdot {}^t Z$ est la matrice A décalée de 1 colonne vers la droite, si bien que l'on retrouve bien la définition précédente pour matrice de Toeplitz. Le rang de déplacement de A est comme auparavant l'entier $\alpha_+(A) = \text{rang}(\phi_+(A))$. On appelle également générateur de déplacement pour l'opérateur $\phi_+$ un couple (G, H) de matrices $n \times \alpha$ vérifiant l'égalité $\phi_+(A) = G \cdot {}^t H$. Si $\alpha_+(A) \ll n$, on dit que A est quasi-Toeplitz. L'entier $\alpha_+(A)$ est appelé longueur du générateur (G, H).

<figure><img src="image_29.png" /></figure>

L'indice + est introduit pour faire pendant à un indice – utilisé ultérieurement (cf Définition 10.8).

Intuitivement, le rang de déplacement mesure à quel point la matrice A est loin d'être de Toeplitz ou d'une concaténation de Toeplitz.

**Exemple 10.4** Les matrices de Toeplitz et Sylvester sont quasi-Toeplitz, de rang de déplacement 2.

Une autre extension permet de définir des matrices *quasi-Vandermonde* et *quasi-Cauchy*, mais pour des opérateurs de déplacement différents :

— $V_a \mapsto V_a - \text{Diag}_a \cdot V_a \cdot {}^t Z$ pour le cas Vandermonde;

— $C_{a,b} \mapsto C_{a,b} - \text{Diag}_a^{-1} \cdot C_{a,b} \cdot \text{Diag}_b$ pour le cas Cauchy.

Ici, pour deux vecteurs $\mathbf{a} = (a_0, \dots, a_{n-1})$ et $\mathbf{b} = (b_0, \dots, b_{n-1})$, on désigne par $V_a$ la matrice de Vandermonde $(a_i^j)_{i,j=0}^{n-1}$, par $\text{Diag}_a$ la matrice diagonale, dont la diagonale porte le vecteur $\mathbf{a}$, et par $C_{a,b}$ la matrice de Cauchy $(1/(a_i - b_j))_{i,j=0}^{n-1}$.

Dans tous ces cas, l'application $\phi_{M,N}$ est définie par $\phi_{M,N}(A) = A - \text{MAN}$, avec M et N bien choisies (en fonction de la structure visée), et le rang de déplacement de A est défini comme étant le rang de la matrice $\phi_{M,N}(A)$.

**Exercice 10.2** Calculer les rangs de déplacement des matrices de Vandermonde et de Cauchy, pour les opérateurs de déplacement définis ci-dessus. ■
Résultat principal

L'idée-clé derrière les algorithmes rapides pour les matrices structurées de rang de déplacement $\alpha$ est l'utilisation des générateurs de déplacement comme structure de données compacte, de taille $O(\alpha n)$, donc proportionnelle à la taille $n$ de la matrice, dans le cas où $\alpha \ll n$. Le résultat principal de ce chapitre est le suivant.

**Théorème 10.2** (Cadre quasi-Toeplitz, quasi-Vandermonde, quasi-Cauchy). Soit $\alpha$ le rang de déplacement adéquat, A une matrice de $M_n(\mathbb{K})$ donnée par des générateurs de déplacement de taille $n \times \alpha$, et soit $b$ un vecteur de $\mathbb{K}^n$. Alors, il est possible de :

1. calculer le déterminant de A;

2. calculer le rang de A;

3. calculer une solution du système $Ax = b$, ou prouver qu'il n'y en a aucune en $\mathbb{O}(\alpha^2 n)$ opérations dans $\mathbb{K}$. Plus exactement, cette complexité s'exprime en termes de la fonction de multiplication polynomiale, et vaut :

1. $O(\alpha^2 M(n) \log n)$ dans le cas quasi-Toeplitz;

2. $O(\alpha^2 M(n) \log^2 n)$ dans les cas quasi-Vandermonde et quasi-Cauchy.

Ce théorème a pour conséquence une *algorithmique unifiée* pour résoudre en complexité quasi-optimale différents problèmes sur les polynômes et les séries.

**Corollaire 10.3** On peut calculer en $O(M(n) \log n)$ opérations arithmétiques :

1. le pgcd étendu et le résultant de deux polynômes de degré borné par $n$;

2. un approximant de Padé de type $(n, n)$ d'une série tronquée à précision $2n$.

*Démonstration.* Le résultant de deux polynômes $A, B \in \mathbb{K}[X]$ s'obtient comme le déterminant de leur matrice de Sylvester $Syl(A, B)$, qui a un rang de déplacement au plus 2. Le degré du pgcd $G = \text{pgcd}(A, B)$ s'obtient par un calcul de rang de la matrice de Sylvester : $\text{deg}(G) = \text{deg}(A) + \text{deg}(B) - \text{rang}(\text{Syl}(A, B))$. Une fois connu le degré du pgcd, la relation de Bézout $UA + VB = G$, avec les contraintes $\text{deg}(U) < \text{deg}(B) - \text{deg}(G)$ et $\text{deg}(V) < \text{deg}(A) - \text{deg}(G)$, se traduit en un système linéaire en les coefficients de $U, V$ et $G$, dont la matrice est une concatenation de trois matrices de Toeplitz, donc une matrice quasi-Toeplitz de rang de déplacement au plus 3. Des considérations similaires s'appliquent pour le calcul d'approximants de Padé. $\blacksquare$

**Corollaire 10.4** Étant données $n$ series $f_1, \dots, f_n$ de $\mathbb{K}[[X]]$ connues à précision $\sigma = \sum_i (d_i + 1) - 1$, il est possible d'en calculer un approximant de Padé-Hermite $(p_1, \dots, p_n)$ de type $(d_1, \dots, d_n)$ en $O(n^2 M(\sigma) \log \sigma)$ opérations dans $\mathbb{K}$.

*Démonstration.* Il s'agit d'un problème linéaire en les coefficients de l'approximant cherché. La matrice est quasi-Toeplitz, de rang de déplacement borné par $n$. $\blacksquare$

Dans la suite, nous allons esquisser les grandes lignes de la preuve du Théorème 10.2 dans le cas particulier où la matrice A est inversible et « suffisamment générique » (tous les mineurs non nuls), et uniquement dans le cas quasi-Toeplitz.
## 10.2 Le cas quasi-Toeplitz

Dans le reste du chapitre, nous allons nous concentrer uniquement sur le cas quasi-Toeplitz, car il contient les principales idées algorithmiques, et il couvre beaucoup d'applications ; les cas quasi-Vandermonde et quasi-Cauchy s'y ramènent et sont plus techniques.

Nous allons montrer que la notion de matrice quasi-Toeplitz est une bonne notion de structure, car elle répond aux propriétés suivantes :

(P1) on peut effectuer quasi-optimalement le produit d'une matrice quasi-Toeplitz par un vecteur;

(P2) la somme et le produit de deux matrices quasi-Toeplitz restent quasi-Toeplitz;

(P3) l'inverse aussi.

Or, ces propriétés forment le minimum nécessaire pour concevoir un algorithme de type inversion de Strassen (page 173) dans la représentation compacte par générateurs de déplacement.

### Produit matrice-vecteur, cas quasi-Toeplitz

Pour montrer la propriété (P1) ci-dessus, le point clé est le suivant.

**Proposition 10.5** — **Formule ΣLU.** L'opérateur $\phi_+$ : $A \mapsto A-Z \cdot A \cdot {}^tZ$ est inversible. Plus exactement, on a la formule suivante, appelée *formule ΣLU* :

$$A - Z \cdot A \cdot {}^tZ == G \cdot {}^tH \quad \text{si et seulement si} \quad A = \sum_{i=1}^{\alpha} L(x_i) \cdot U(y_i),$$

où les $x_i$ (resp. $y_i$) sont les colonnes du générateur G (resp. H), et où, pour un vecteur colonne $v = {}^t[v_0 \cdots v_{n-1}]$, on note $L(v)$ la matrice de Toeplitz triangulaire inférieure

$$\begin{pmatrix} v_0 & 0 & \cdots & 0 \\ v_1 & v_0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ v_{n-1} & \cdots & v_1 & v_0 \end{pmatrix}$$

et $U(v)$ la matrice de Toeplitz triangulaire supérieure $^tL(v)$.

*Démonstration.* Par linéarité, il suffit de traiter le cas $\alpha = 1$. Si $C = L(a)U(b)$, alors un calcul immédiat montre que $c_{i,j} = a_i b_j + a_{i-1} b_{j-1} + \cdots$, et donc $c_{i,j} - c_{i-1,j-1} = a_i b_j$ et $\phi_+(C) = (a_i b_j)^{n-1}_{i,j=0} = a \cdot {}^t b$.

L'implication inverse découle de l'injectivité de $\phi_+$ et est laissée en exercice. ■

La Proposition 10.5 permet de donner une définition équivalente pour le rang de déplacement.

**Définition 10.7** Le nombre $\alpha_+(A)$ est le plus petit entier $\alpha$ tel qu'il existe une
décomposition de la forme

$$A = \sum_{i=1}^{\alpha} L_i U_i$$

pour des matrices $L_i$ de Toeplitz inférieures, et $U_i$ de Toeplitz supérieures.

**Exemple 10.5** Si $A$ est de Toeplitz, alors en notant $A_{\text{inf}}$ sa partie inférieure et $A_{\text{ssup}}$ sa partie strictement supérieure, on a que $A$ admet la décomposition $A = A_{\text{inf}} \cdot I_n + I_n \cdot A_{\text{ssup}}$, donc $\alpha_+(A) \le 2$

La Définition 10.7 permet de prouver la propriété (P1) en page 201.

**Corollaire 10.6** Si $A$ est donnée en représentation compacte par une paire de générateurs de déplacement (G, H) de taille $n \times \alpha$, alors le produit matrice-vecteur $Av$ peut s'effectuer en $O(\alpha M(n))$ opérations arithmétiques.

*Démonstration.* Le produit $Av$ s'écrit comme la somme de $\alpha$ termes, tous de la forme $L(x_i)(U(y_i)v)$. Or, chacun de ces termes peut être calculé en $O(M(n))$ opérations, grâce au Lemme 10.1. ■

## Addition et produit en représentation compacte par générateurs

Pour justifier la propriété (P2) en page 201, on va montrer plus, à savoir que la représentation par générateurs permet un calcul efficace (en temps quasi-linéaire) des opérations d'additions et de multiplication de deux matrices quasi-Toeplitz.

**Proposition 10.7** — Opérations matricielles en représentation compacte.

Soient (T, U) générateur de déplacement de longueur $\alpha$ pour A, et (G, H) générateur de longueur $\beta$ pour B. Alors

1. $([T|G],[U|H])$ sont des générateurs de déplacement pour $A+B$ de longueur $\alpha+\beta$;

2. $([T|W|a],[V|H|-b])$ sont des générateurs de déplacement pour AB, de longueur $\alpha+\beta+1$,

où $V := {}^tB \cdot U$, $W := Z \cdot A \cdot {}^tZ \cdot G$, et où le vecteur $a$ (resp. $b$) est la dernière colonne de $Z \cdot A$ (resp. de $Z \cdot {}^tB$).

La preuve, basée sur une vérification immédiate, est laissée en exercice.

**Corollaire 10.8** En représentation compacte par générateurs de déplacement, de longueurs au plus $\alpha$ pour A et B, on peut calculer

1. la somme $A+B$ en $O(\alpha n)$ opérations dans $\mathbb{K}$;

2. le produit $AB$ en $\text{Mul}(n, \alpha) := O(\alpha^2 M(n))$ opérations dans $\mathbb{K}$.

*Démonstration.* Le seul point non trivial est le calcul des matrices $V$ et $W$ et des vecteurs $a$ et $b$. Tout repose sur la formule $\Sigma LU$. Si $B$ s'écrit $\sum_{i=1}^{\beta} L(x_i)U(y_i)$, alors
sa transposée s'écrit $^tB = \sum_{i=1}^{\beta} L(y_i)U(x_i)$, et le calcul de $V = {}^tB \cdot U$ se ramène à $2\alpha\beta$ multiplications polynomiales, chacune pouvant s'effectuer en $O(M(n))$. Le même raisonnement s'applique au calcul de $W$. Enfin, le calcul de $a$ revient à multiplier $A$ par la colonne $^t[0, \dots, 0, 1]$ et une observation similaire pour $b$ permet de conclure. ■

### Inversion en représentation compacte par générateurs

Afin de prouver la propriété (P3) en page 201, il est commode d'introduire un opérateur dual de déplacement.

**Définition 10.8** L'opérateur de $\phi_-$-déplacement d'une matrice A de $M_n(\mathbb{K})$ est défini par l'égalité

$$ \begin{aligned} \phi_-(A) &= A - {}^tZ \cdot A \cdot Z \\ &= A - (A \text{ décalée de 1 vers le haut et vers la gauche}). \end{aligned} $$

Le rang de déplacement $\alpha_-$ est défini par la formule

$$ \alpha_-(A) = \mathrm{rang}(\phi_-(A)). $$

On peut montrer (de la même façon que pour $\phi_+$) que l'on dispose d'une caractérisation équivalente pour $\phi_-$, et que $\phi_+$ et $\phi_-$ sont liés.

### Lemme 10.9

1. $\alpha_-(A)$ est le plus petit entier $\alpha$ tel que $A$ puisse s'écrire sous la forme $\sum_{i=1}^{\alpha} U_i L_i$, avec les $L_i$ de Toeplitz inférieures, et les $U_i$ de Toeplitz supérieures.

2. On a la formule $\Sigma UL$ :

$$ A - {}^tZ \cdot A \cdot Z = \sum_{i=1}^{\alpha} x_i \cdot {}^t y_i \quad \text{si et seulement si} \quad A = \sum_{i=1}^{\alpha} U(\mathrm{rev}(x_i)) \cdot L(\mathrm{rev}(y_i)). $$

Ici, pour $v = {}^t[v_0, \dots, v_{n-1}]$, on note $\mathrm{rev}(v) = {}^t[v_{n-1}, \dots, v_0]$.

3. Une matrice de Toeplitz pour $\phi_+$ est une matrice de Toeplitz pour $\phi_-$.

La preuve de ce résultat découle de la proposition suivante.

**Proposition 10.10** — Conversion $\Sigma LU \leftrightarrow \Sigma UL$. Pour toute matrice A, l'inégalité $|\alpha_+(A) - \alpha_-(A)| \le 2$ est satisfaite. De plus, on peut effectuer les conversions d'une représentation $\Sigma LU$ à une représentation $\Sigma UL$, et inversement, en $O(\alpha M(n))$ opérations dans $\mathbb{K}$.

*Démonstration.* Il suffit de prouver l'identité

$$ L(x) \cdot U(y) = I_n \cdot L(y') + U(x') \cdot I_n - U(x'') \cdot L(y''), $$

où $x'' = Z \cdot \mathrm{rev}(x)$, $y'' = Z \cdot \mathrm{rev}(y)$, $y' = \mathrm{rev}( {}^tU(y) \cdot {}^tL(x) \cdot f)$, $x' = \mathrm{rev}(L(x) \cdot U(y) \cdot f)$, et $f = {}^t[0, \dots, 0, 1]$. ■
**Entrée** Une matrice « générique » $A \in M_n(\mathbb{K})$, de taille $n = 2^k$, représentée par des générateurs (pour l'opérateur de déplacement $\phi_+$).

**Sortie** Son inverse $A^{-1}$, représentée par des générateurs (pour l'opérateur de déplacement $\phi_-$).

1. Si $n = 1$, renvoyer $A^{-1}$.

2. Calculer des $\phi_+$-générateurs pour $a, b, c, d \in M_{n/2}(\mathbb{K})$,

$$ \text{où } A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}. $$

3. Calculer récursivement des $\phi_-$-générateurs pour $e := a^{-1}$.

4. Calculer des $\phi_+$-générateurs pour $S := d - ceb$.

5. Calculer récursivement des $\phi_-$-générateurs pour $t := S^{-1}$.

6. Renvoyer des $\phi_-$-générateurs de $A^{-1} = \begin{pmatrix} x & y \\ z & t \end{pmatrix}$ via les formules de Strassen $y := -ebt, z := -tce$ et $x := e + ebtce$.

**Algorithme 10.1** – Algorithme de type Strassen pour inverser une matrice quasi-Toeplitz.

**Théorème 10.11** Soit $A \in M_n(\mathbb{K})$ une matrice quasi-Toeplitz inversible. Alors son inverse est aussi quasi-Toeplitz et $\alpha_+(A^{-1}) = \alpha_-(A)$.

**Démonstration.** On a la suite d'égalités

$$ \begin{align*} \alpha_-(A) &= \mathrm{rang}(A - {}^tZ \cdot A \cdot Z) = \mathrm{rang}(I_n - A^{-1} \cdot {}^tZ \cdot A \cdot Z) = \mathrm{rang}(I_n - Z \cdot A^{-1} \cdot {}^tZ \cdot A) \\ &= \mathrm{rang}(A^{-1} - Z \cdot A^{-1} \cdot {}^tZ) = \alpha_+(A^{-1}), \end{align*} $$

dans laquelle nous avons utilisé un fait d'algèbre linéaire élémentaire : si $A, B$ sont deux matrices quelconques, alors $\mathrm{rang}(I_n - A \cdot B) = \mathrm{rang}(I_n - B \cdot A)$. ■

## Résolution rapide de systèmes quasi-Toeplitz

Le Théorème 10.11 ne nous dit pas *comment* calculer des générateurs pour $A^{-1}$. Cela est fait dans le dernier résultat de ce chapitre.

**Théorème 10.12** Soit $A \in M_n(\mathbb{K})$ une matrice quasi-Toeplitz inversible, de rang de déplacement $\alpha$, donnée par des générateurs $G$ et $H$ de longueur $\alpha$. On peut calculer des générateurs de longueur $\alpha$ pour $A^{-1}$ en $O(\alpha^2 M(n) \log n)$ opérations dans $\mathbb{K}$.

À partir de cette représentation de l'inverse et d'un vecteur $b \in \mathbb{K}^n$, le système $Ax = b$ peut se résoudre pour un supplément de $O(\alpha M(n))$ opérations dans $\mathbb{K}$.

**Démonstration.** L'idée est d'adapter au cas structuré l'algorithme d'inversion de Strassen présenté au Chapitre 8 (page 173). La différence réside dans le choix de la structure
de données utilisée pour représenter et calculer avec les matrices quasi-Toeplitz. Ici, on utilise comme structure de données compacte les générateurs de déplacement pour les opérateurs $\phi_+$ et $\phi_-$ (ou, de façon équivalente, les représentations $\Sigma$LU et $\Sigma$UL).

La méthode, de type « diviser pour régner », est résumée dans l'Algorithme 10.1. Sa correction est héritée de celle de l'algorithme d'inversion de Strassen ; l'hypothèse sur la généricité de $A$ est utilisée pour garantir que tous les mineurs qu'on doit inverser au cours de l'algorithme sont bien inversibles.

Sa complexité arithmétique $C(n)$ obéit à la récurrence $C(1) = 1$ et, si $n \ge 2$,

$$C(n) = 2C(n/2) + O(\mathrm{Mul}(n, \alpha)) + O(\alpha^2 n + \alpha M(n)),$$

où le coût en $O(\alpha^2 n + \alpha M(n))$ provient des conversions entre les représentations $\phi_+$ et $\phi_-$ et des calculs de minimisation de longueurs de générateurs, et la notation $\mathrm{Mul}(n, \alpha)$ désigne la complexité du produit matriciel pour des matrices $n \times n$ données par générateurs de taille $n \times \alpha$. La preuve se conclut en utilisant l'estimation

$$\mathrm{Mul}(n, \alpha) = O(\alpha^2 M(n)),$$

$$\text{qui entraîne } C(n) = O(\alpha^2 M(n) \log n). \quad \blacksquare$$

Exercices

**Exercice 10.3** Soit $A \in \mathcal{M}_n(\mathbb{K})$ une matrice quasi-Toeplitz de rang de déplacement $\alpha \ll n$, représentée de façon compacte par des générateurs $(G, H)$ de taille $n \times \alpha$, c'est-à-dire telle que $\phi_+(A) = G \cdot {}^t H$.

1. Soient $v_1, \dots, v_\alpha$ des vecteurs quelconques de $\mathbb{K}^n$. Montrer que le calcul de tous les produits $A \cdot v_\ell$, $1 \le \ell \le \alpha$, se ramène, grâce à la formule $\Sigma$LU, au problème suivant :

(P) Étant donnés des polynômes $G_j, H_j, V_j \in \mathbb{K}[X]$ ($j \le \alpha$) de degrés au plus $n-1$,

$$\text{calculer } A_{\ell} = \sum_{j=1}^{\alpha} G_j (H_j V_{\ell} \bmod X^n), \quad 1 \le \ell \le \alpha.$$

2. Donner un premier algorithme qui résout le problème (P) et estimer sa complexité.

3. Montrer que le problème (P) admet la reformulation matricielle suivante :
   (MP) Étant données des matrices polynomiales $G, V$ dans $\mathcal{M}_{\alpha \times 1}(\mathbb{K}[X])$ et $H$ dans $\mathcal{M}_{1 \times \alpha}(\mathbb{K}[X])$, toutes de degrés au plus $n-1$, calculer $(VH \bmod X^n)G$.

4. Soient $A, B$ et $C$ des matrices polynomiales de tailles $(n \times p)$, $(p \times n)$ et $(n \times p)$ et de degré au plus $d$. Montrer que le produit $ABC$ peut être calculé en $O(\frac{n}{p} \mathrm{MM}(p, d))$ opérations dans $\mathbb{K}$.

5. Proposer un algorithme de type « diviser pour régner » (par rapport à $n$) résolvant le problème (MP) en complexité $O(\frac{1}{\alpha} \mathrm{MM}(\alpha, n))$.

6. Conclure qu'on peut multiplier, dans la représentation par générateurs de déplacement, deux matrices quasi-Toeplitz de $\mathcal{M}_n(\mathbb{K})$ de rang de déplacement au plus $\alpha$ en $O(\frac{1}{\alpha} \mathrm{MM}(\alpha, n))$ opérations dans $\mathbb{K}$.

■
## Notes

Levinson [Lev47] a donné le premier algorithme de complexité quadratique pour la résolution d'un système linéaire dont la matrice est de Toeplitz, symétrique et définie positive. Trench [Tre64] a par la suite montré que tout l'inverse d'une telle matrice peut se calculer en la même complexité. Ces deux algorithmes sont décrits dans le livre de Golub et Van Loan [GV13, §4.7].

L'algorithme de Trench [Tre64] repose sur une formule, reprise et améliorée par Gohberg et Semencul [GS72], qui fournit une représentation explicite ΣLU à deux termes de l'inverse d'une matrice de Toeplitz A, les sommands étant construits à partir de la première ligne et de la première colonne de $A^{-1}$. Brent, Gustavson et Yun [BGY80; GY79] ont montré comment ramener le calcul des premières ligne et colonne de $A^{-1}$ à un calcul d'approximation de Padé. Ils ont aussi mis en évidence les liens profonds entre la résolution des systèmes de Toeplitz, l'algorithme d'Euclide étendu et l'approximation de Padé, et ont proposé le premier algorithme pour la résolution d'un système défini par une matrice de Toeplitz inversible quelconque de taille *n* en $O(M(n) \log n)$ opérations arithmétiques.

Parallèlement, la notion de rang de déplacement a été introduite par Kailath, Kung et Morf dans [KKM79b], pour pouvoir résoudre un système quasi-Toeplitz inversible de rang de déplacement $\alpha$ en coût $O(\alpha n^2)$ [KKM79a]. Des versions rapides, de complexité $O(\alpha^d M(n) \log n)$ ont été obtenues indépendamment par Bitmead et Anderson [BA80] (avec $d = 4$) et Morf [Mor80] (avec $d = 2$), sous des hypothèses de forte régularité sur l'entrée. Ces hypothèses ont été ensuite levées par Kaltofen [Kal94; Kal95b]. Les extensions au cas quasi-Cauchy sont dues à Cardinal [Car99], Pan et Zheng [PZ00]. Le cas quasi-Vandermonde est traité dans des articles de Pan [Pan90], et de Gohberg et Olshevsky [GO94]. La meilleure complexité asymptotique vis-à-vis simultanément de la taille de la matrice et de son rang de déplacement est due à Bostan, Jeannerod et Schost [BJS08], qui ont trouvé comment introduire du produit rapide de matrices denses dans l'algorithmique des matrices structurées, et ont réduit le coût à $\tilde{O}(\alpha^{\theta-1}n)$ pour la résolution des systèmes quasi-Toeplitz, quasi-Vandermonde et quasi-Cauchy.

Une bonne référence sur les matrices structurées est le livre de Pan [Pan01].

## Bibliographie

* BA80 BITMEAD, Robert R. et Brian D. O. ANDERSON (1980). « Asymptotically fast solution of Toeplitz and related systems of linear equations ». In: *Linear Algebra and its Applications*, vol. 34, p. 103-116.

* BGY80 BRENT, Richard P., Fred G. GUSTAVSON et David Y. Y. YUN (1980). « Fast solution of Toeplitz systems of equations and computation of Padé approximants ». In: *Journal of Algorithms*, vol. 1, n°3, p. 259-295.

* BJS08 BOSTAN, Alin, Claude-Pierre JEANNEROD et Éric SCHOST (2008). « Solving structured linear systems with large displacement rank ». In: *Theoretical Computer Science*, vol. 407, n°1-3, p. 155-181.

* Car99 CARDINAL, Jean-Paul (1999). « On a property of Cauchy-like matrices ». In: *Comptes Rendus de l'Académie des Sciences. Série I. Mathématique*, vol. 328, n°11, p. 1089-1093.
GO94 GOHBERG, I. et V. OLSHEVSKY (1994). « Complexity of multiplication with vectors for structured matrices ». In : *Linear Algebra and its Applications*, vol. 202, p. 163–192.

GS72 GOHBERG, I. C. et A. A. SEMENCUL (1972). « On the inversion of finite Toeplitz matrices and their continuous analogues (in Russian) ». In : *Matematicheskie Issledovaniya*, vol. 7, n°2, p. 201–223.

GV13 GOLUB, Gene H. et Charles F. VAN LOAN (2013). *Matrix computations*. Fourth. Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD.

GY79 GUSTAVSON, Fred G. et David Y. Y. YUN (1979). « Fast algorithms for rational Hermite approximation and solution of Toeplitz systems ». In : *IEEE Transactions on Circuits and Systems*, vol. 26, n°9, p. 750–755.

Kal94 KALTOFEN, Erich (1994). « Asymptotically fast solution of Toeplitz-like singular linear systems ». In : *ISSAC'94 : International Symposium on Symbolic and Algebraic Computation*. Oxford, United Kingdom : ACM Press, p. 297–304.

Kal95b — (1995). « Analysis of Coppersmith's block Wiedemann algorithm for the parallel solution of sparse linear systems ». In : *Mathematics of Computation*, vol. 64, n°210, p. 777–806.

KKM79a KAILATH, T., S. Y. KUNG et M. MORF (1979). « Displacement ranks of a matrix ». In : *American Mathematical Society. Bulletin. New Series*, vol. 1, n°5, p. 769–773.

KKM79b KAILATH, Thomas, Sun Yuan KUNG et Martin MORF (1979). « Displacement ranks of matrices and linear equations ». In : *Journal of Mathematical Analysis and Applications*, vol. 68, n°2, p. 395–407.

Lev47 LEVINSON, Norman (1947). « The Wiener RMS (root mean square) error criterion in filter design and prediction ». In : *Journal of Mathematics and Physics*, vol. 25, p. 261–278.

Mor80 MORF, M. (1980). « Doubling algorithms for Toeplitz and related equations ». In : *IEEE Conference on Acoustics, Speech, and Signal Processing*, p. 954–959.

Pan01 PAN, Victor Y. (2001). *Structured matrices and polynomials*. Unified superfast algorithms. Birkhäuser Boston Inc.

Pan90 PAN, Victor (1990). « On computations with dense structured matrices ». In : *Mathematics of Computation*, vol. 55, n°191, p. 179–190.

PZ00 PAN, Victor Y. et Ailong ZHENG (2000). « Superfast algorithms for Cauchy-like matrix computations and extensions ». In : *Linear Algebra and its Applications*, vol. 310, n°1-3, p. 83–108.

Tre64 TRENCH, William F. (1964). « An algorithm for the inversion of finite Toeplitz matrices ». In : *Journal of the Society for Industrial and Applied Mathematics*, vol. 12, p. 515–522.
<figure><img src="image_30.png" /></figure>
# 11. Solutions rationnelles de systèmes linéaires à coefficients polynomiaux

## Résumé

L'algorithmique des systèmes linéaires à coefficients des polynômes en une variable est très similaire à celle des fractions rationnelles. En outre, il est important de tirer parti du produit rapide de matrices.

On considère le système linéaire

$$A(X)Y(X) = B(X), \qquad (11.1)$$

où A et B sont donnés et Y est inconnu. A est une matrice $n×n$ de polynômes, régulière (de déterminant non nul) et B est un vecteur de polynômes. De manière équivalente, on peut voir A (ou B) comme un polynôme à coefficients des matrices (ou des vecteurs). Pour fixer les notations, on suppose $\deg A \le d$ et $\deg B < d$.

On notera $M_{m,k}(R)$ l'ensemble des matrices de taille $m×k$ à coefficients dans $R$. On notera $K[X]_d$ l'ensemble des polynômes de degré au plus $d$ à coefficients dans le corps $K$. La fonction $M$ est telle que la multiplication dans $K[X]_d$ coûte au plus $M(d)$ opérations dans $K$. L'exposant $\theta$ est tel que la multiplication de deux matrices $n×n$ à coefficients dans $K$ coûte $MM(n) = O(n^\theta)$ opérations dans $K$.

Nous aurons aussi besoin de produits de matrices de $M_n(K[X]_d)$. Dans le cas le plus général, ce produit est connu pour pouvoir être exécuté en $MM(n, d) = O(n^\theta M(d))$ opérations dans $K$, borne que nous utiliserons dans les estimations de complexité. Lorsque le corps $K$ contient suffisamment de points, ce coût descend par évaluation-interpolation à $O(n^\theta d + n^2 M(d) \log d)$; dans les mêmes conditions, en choisissant des points en progression géométrique, cette complexité peut être encore abaissée à $O(n^\theta d + n^2 M(d))$, voir page 106.
**Entrée** A ∈ M_n(ℝ[X]_d), B ∈ M_{n,1}(ℝ[X]_{d-1}).

**Sortie** le vecteur de fractions rationnelles Y tel que AY = B.

1. Calculer le développement en série de A⁻¹B à précision 2nd.

2. Reconstruire les coefficients de Y par approximant de Padé.

Algorithme 11.1 – Résolution de A(X)Y(X) = B(X).

## 11.1 Des séries aux solutions rationnelles

Une première observation est que la structure de la solution cherchée est donnée par les formules de Cramer.

**Lemme 11.1** Le système possède une solution dont les coordonnées sont des fractions rationnelles. Ses numérateurs et dénominateurs ont degrés bornés par nd - 1 et nd.

*Démonstration.* Le système (11.1) se récrit $A_1y_1 + \cdots + A_ny_n = B$, où $y_i$ est la ième coordonnée de Y et $A_i$ la ième colonne de A. La formule de Cramer

$$ \det(A_1, \dots, A_{i-1}, B, A_{i+1}, \dots, A_n) = y_i \det(A), $$

est obtenue en remplaçant B par sa valeur dans le membre gauche et en développant le déterminant.

Donc $y_i$ est le quotient de déterminants de matrices appartenant à $M_n(ℝ[X]_d)$, qui ont nécessairement degré au plus $nd - 1$ pour le numérateur et $nd$ pour le dénominateur. ■

L'Algorithme 11.1 de résolution est justifié par la complexité quasi-optimale des approximants de Padé (Chapitre 7). Dans tous les algorithmes qui vont suivre, c'est la recherche de série solution qui va dominer la complexité.

## 11.2 Développement comme une fraction rationnelle

La Proposition 3.5 (page 64) montre que la méthode de Newton fonctionne pour toute matrice inversible de séries. Avec cet algorithme, le calcul de la fraction rationnelle solution de (11.1) demande $O(n^θ M(nd))$ opérations dans $ℝ$.

Il est possible d'améliorer l'efficacité lorsque la matrice est une matrice de polynômes. La base de cette amélioration est contenue dans le lemme suivant.

**Lemme 11.2** Soit A(X) ∈ M_n(ℝ[X]_d) et B(X) ∈ M_{n,m}(ℝ[X]_{d-1}), avec A inversible. Pour tout k, il existe une matrice B_k ∈ M_{n,m}(ℝ[X]_{d-1}) telle que

$$ A^{-1}B = a_0 + a_1X + \cdots + a_{k-1}X^{k-1} + X^k A^{-1}B_k, \quad (11.2) $$

où $a_i \in M_{n,m}(ℝ)$.
**Entrée** A, B, k et S = A⁻¹ mod Xᵏ.

**Sortie** a₀, ..., aₖ₋₁ et Bₖ définis par l'équation (11.2).

1. Calculer SB =: a₀ + ... + aₖ₋₁Xᵏ₋₁ mod Xᵏ.

2. Calculer Bₖ = (B - A(a₀ + ... + aₖ₋₁Xᵏ₋₁))X⁻ᵏ.

3. Renvoyer a₀, ..., aₖ₋₁, Bₖ.

## Algorithme 11.2 – Développement de A⁻¹B.

**Démonstration.** Si les aᵢ sont les coefficients du développement en série de A⁻¹B, alors

$$B - A(a_0 + \cdots + a_{k-1}X^{k-1})$$

est une matrice de Mₙ,ₘ(K[X]ᵢⱼ₋₁) qui, par construction, est divisible par Xᵏ, d'où le lemme. ■

Ce résultat se traduit par l'Algorithmme 11.2. La proposition suivante estime le coût de cet algorithme. Elle peut être vue comme l'analogue matriciel du Théorème 4.9 (page 88).

**Proposition 11.3** Soit A ∈ Mₙ(K[X]ᵢⱼ) avec A(0) inversible, alors on peut calculer les N ≥ d premiers coefficients de A⁻¹ en O(n°N M(d)/d) opérations dans K.

Pour B ∈ Mₙ,₁(K[X]ᵢⱼ₋₁), on peut calculer les N ≥ d premiers coefficients de A⁻¹B en O(n²N M(d)/d) opérations dans K.

**Démonstration.** L'algorithme calcule d'abord l'inverse S = A⁻¹ mod Xᵈ par l'algorithme de Newton, en O(MM(n, d)) = O(n° M(d)) opérations dans K.

Ensuite, on applique N/d fois l'algorithme ci-dessus avec k = d pour calculer à chaque itération d coefficients et un nouveau Bᵢⱼ₋₁. Si on part de B₀ = I, le résultat fournit les N premiers coefficients de A⁻¹; si B₀ = B, alors on obtient les coefficients de A⁻¹B.

Les deux étapes de l'algorithme (avec k = d) coûtent MM(n, d) = O(n° M(d)) opérations dans K si B a n colonnes, O(n² M(d)) s'il n'en a qu'une. ■

Avec cet algorithme, le calcul de la fraction rationnelle solution de (11.1) demande O(n³ M(d)) opérations dans K.

## 11.3 L'algorithme de Storjohann

L'algorithme de Storjohann permet de calculer les N premiers coefficients du développement en série de A⁻¹B en O(n°⁻¹N log N M(d)) opérations dans K. Si θ < 3, cette quantité croît avec n moins vite que la taille de l'inverse A⁻¹, qui n'est donc pas calculée en entier. Ainsi, la résolution du système linéaire (11.1) a un coût en O(n° M(d) log(nd)) opérations dans K.
**Entrée** A, B, k et S = $[A^{-1}]_{k-2d+1}^{2d-2}$.

**Sortie** B_k défini par l'équation (11.2).

1. Calculer U := [SB]_{k-d}^{d-1}.

2. Calculer $B_k := [B - AU]_k^{d-1}X^{-k}$.

3. Renvoyer $B_k$.

Algorithme 11.3 – Développement de $A^{-1}B$ tronqué.

L'algorithme repose sur le développement de $A^{-1}B$ de la section précédente, joint d'une part à une technique de type « diviser pour régner », d'autre part au regroupe- ment des calculs intermédiaires pour remplacer des groupes de $n$ produits matrice- vecteur par des produits matrice-matrice.

Pour commencer, on peut modifier l'algorithme de développement de $A^{-1}B$ de la section précédente pour calculer $B_k$ sans calculer tous les coefficients $a_0, \dots, a_{k-1}$. L'entrée nécessaire est moins grosse, et la complexité plus faible lorsque $k$ est grand devant $d$. Pour une série $V = v_0 + v_1 X + \dots$, on note

$$[V]_a^b := v_a X^a + \dots + v_{a+b} X^{a+b},$$

avec la convention que les coefficients d'indice négatif de $V$ sont nuls. Cela donne l'Algorithme 11.3.

*Démonstration.* Pour prouver la correction de cet algorithme, il faut s'assurer que les troncatures de séries sont suffisantes pour calculer le même polynôme $B_k$ que précédemment.

Pour la première étape de l'algorithme, il suffit d'observer que $B$ ayant degré au plus $d-1$, le coefficient de $X^i$ dans $A^{-1}B$ ne dépend que de $[A^{-1}]_{i-d-1}^{d+1}$, pour tout $i$. Donc les coefficients calculés par cette première étape sont les mêmes que les $a_i$ que précédemment, pour $i = k-d, \dots, k-1$.

Ensuite, il faut calculer $[X^k B_k]_k^{d-1}$. L'extraction des coefficients de l'équation (11.2) donne

$$[X^k B_k]_k^{d-1} = [B - A(a_0 + \dots + a_{k-1} X^{k-1})]_k^{d-1} = [B - A(a_{k-d} X^{k-d} + \dots + a_{k-1} X^{k-1})]_k^{d-1},$$

ce qui conclut la preuve.

<figure><img src="image_31.png" /></figure>

Une observation importante concernant cet algorithme est que c'est le même polynôme $S$ qui peut servir pour calculer $B_{i+k}$ pour tout $i$. L'idée est alors de grouper plusieurs vecteurs $B_i$ et de calculer les $B_{i+k}$ correspondants par des produits matrice-matrice. Noter que la ressemblance entre cette idée et celle de l'algorithme de Keller-Gehrig en page 176.

Ainsi, si l'on connaît $B_0, B_{2m}, \dots, B_{2sm}$ et $[A^{-1}]_{m-2d}^{2d}$, alors le calcul des vecteurs suivants $B_m, B_{3m}, \dots, B_{(2s+1)m}$ ne requiert que $O(n^{\theta-1} s M(d))$ opérations dans $\mathbb{K}$.

En itérant cette idée, en partant de $B_0$ et en supposant connus $[A^{-1}]_{2k-2d+1}^{2d-2}$ pour $k=0, \dots, \lfloor \log(N/d) \rfloor := k_{\max}$, on obtient d'abord $B_0, B_{2k_{\max}}$, puis à l'étape suivante $B_0, B_{2k_{\max}-1}, B_{2k_{\max}}$, $B_{3,2k_{\max}-1}$, et ainsi de suite jusqu'à calculer toute la suite $B_0, B_d, B_{2d}, \dots$,
**Entrée** $S = [A^{-1}]_{0}^{d-1}$, $T = [A^{-1}]_{2k-2d+1}^{2d-2}$ et $B_{2k-d}$ défini par (11.2) avec $B = I$.

**Sortie** $B_{2k+1-d}$ et $[A^{-1}]_{2k+1-2d+1}^{2d-2}$.

1. Calculer $[A^{-1}]_{2k+1-2d+1}^{d-2} = [TB_{2k-d}]_{2k+1-2d+1}^{d-2}$.

2. Calculer $B_{2k+1-d} := [-ATB_{2k-d}]_{2k}^{d-1}/X^{2k}$.

3. Calculer $[A^{-1}]_{2k+1-d}^{d-1} = [X^{2k+1-d}B_{2k+1-d}S]_{2k+1-d}^{d-1}$.

4. Renvoyer $B_{2k+1-d}$ et $[A^{-1}]_{2k+1-2d+1}^{2d-2}$.

**Algorithme 11.4** – Développement de $A^{-1}$ — indices puissances de 2.

$B_{d[N/d]}$ en $O(k_{max} n^{\theta-1} N M(d)/d)$ opérations dans $\mathbb{K}$. En multipliant alors ces vecteurs par $[A^{-1}]_{0}^{d}$ on obtient finalement les N premiers coefficients de $A^{-1}B$ pour le même coût.

Il reste à voir comment calculer les $[A^{-1}]_{2k-2d-1}^{2d-2}$. Là encore, le point de départ est l'identité (11.2), avec $B = I$, $k = m$ et $k = p$ :

$$ 
\begin{aligned}
A^{-1} &= a_0 + \cdots + a_{m-1} X^{m-1} + X^m A^{-1} B_m \\
&= a_0 + \cdots + a_{m-1} X^{m-1} + X^m (a_0 + \cdots + a_{p-1} X^{p-1} + X^p A^{-1} B_p) B_m.
\end{aligned}
$$

La seconde ligne est obtenue par substitution de la valeur de $A^{-1}$ donnée par la première ligne avec $m = p$. Ces équations entraînent pour tout $\ell \ge 0$ tel que $m+p-\ell \ge d$

$$ 
\left\{ 
\begin{array}{ll}
B_{m+p} & = [-A[A^{-1}]_{p-2d+1}^{2d-2} B_m]_{p}^{d-1} X^{-p}, \\
[A^{-1}]_{m+p-\ell}^{\ell-1} & = [[A^{-1}]_{p-\ell-d+1}^{\ell+d-2} B_m]_{m+p-\ell}^{\ell-1}
\end{array}
\right.
$$

L'Algorithmme 11.4 s'en déduit en utilisant cette identité avec $m = 2^k - d$, $p = 2^k$ et $\ell = d - 1$. Pour résumer, ces algorithmes mènent au résultat suivant.

**Théorème 11.4** — Storjohann, 2002. Soient A une matrice $n \times n$ polynomiale de degré $d$ avec $A(0)$ inversible et B un vecteur polynomial de degré au plus $d-1$, alors on peut calculer le numérateur et le dénominateur de $A^{-1}B$ en $O(n^\theta M(d) \log(nd))$ opérations dans $\mathbb{K}$.

L'énoncé suivant est un analogue du théorème 11.4 dans le cas des matrices à coefficients entiers. Sa preuve est toutefois bien plus technique et sera omise.

**Théorème 11.5** — Storjohann, 2005. Soit A une matrice $n \times n$ inversible contenant des entiers de taille binaire au plus $b$ et soit B un vecteur $n \times 1$ contenant des entiers de taille $O(b)$. On peut calculer le numérateur et le dénominateur de $A^{-1}B$ en $O(n^\theta \log(n) M(d) \log(d))$ opérations binaires, où $d = b + \log n$.
## Notes

La conception d'algorithmes pour résoudre des systèmes linéaires à coefficients polynomiaux (resp. entiers) par développement de Taylor (resp. p-adique) et reconstruction rationnelle remonte aux travaux de Moenck et Carter [MC79] (dont sont issus l'Algorithme 11.1 et la Proposition 11.3) d'une part, et de Dixon [Dix82] d'autre part.

Dans tout ce texte, nous avons supposé que $A(0)$ est inversible. En fait, les résultats s'étendent au cas où $A$ est inversible sans que $A(0)$ le soit. Il suffit de tirer aléatoirement un point et d'effectuer les développements en série au voisinage de ce point. L'algorithme devient probabiliste. Si la caractéristique n'est pas nulle, il se peut que $A$ soit inversible sans que sa valeur en aucun point du corps ne le soit. On peut alors construire une extension du corps assez grande pour trouver un point où effectuer ces calculs. Ces considérations ont été prises en compte dans des travaux de Storjohann [Sto02; Sto03].

En 2003, Storjohann [Sto03] a montré que le déterminant d'une matrice polynomiale de taille $n$ et degré au plus $d$ pouvait se calculer en $O(n^\theta M(d) \log^2 n)$ opérations arithmétiques. Storjohann et Villard [SV05] sont ensuite parvenus à calculer le rang et une base du noyau d'une telle matrice en $\tilde{O}(n^\theta d)$ opérations arithmétiques.

Pour $n$ une puissance de 2, Jeannerod et Villard [JV05] ont donné un algorithme déterministe qui calcule l'inverse d'une matrice polynomiale inversible de taille $n$ et degré au plus $d$ génériquement en $\tilde{O}(n^3 d)$ opérations arithmétiques. Ils ont également étudié les réductions entre diverses opérations sur les matrices polynomiales en collaboration avec Giorgi [GJV03]. Plus récemment, Zhou, Labahn et Storjohann [ZLS15] ont amélioré le résultat de Jeannerod et Villard [JV05], en levant l'hypothèse de généricité et celle sur la taille $n$.

Storjohann a étendu tous les algorithmes de ce chapitre et tous ses autres résultats de 2002 et 2003 [Sto02; Sto03] au cas entier dans son article de 2005 [Sto05]. Il est le premier à être parvenu à montrer que le déterminant d'une matrice $n \times n$ formée d'entiers de taille binaire au plus $b$ peut se calculer, à des facteurs logarithmiques près, aussi vite que le coût de la multiplication de deux telles matrices, plus exactement en $O(n^\theta \log(n)^2 M(d) \log(d))$ opérations binaires, où $d = b+\log n$. En particulier, l'exposant du calcul du déterminant d'une matrice entière est bien égal à l'exposant $\omega$ du produit matriciel. Le meilleur algorithme précédemment connu, dû à Kaltofen et Villard [KV04] avait comme exposant l'expression exotique $\omega + \frac{1-\zeta}{\omega^2-(2+\zeta)\omega+2}$, où $\zeta \in [0,1]$ est un réel tel que la multiplication de deux matrices de taille $n \times n$ et $n \times n^\zeta$ puisse s'effectuer en $O(n^{2+o(1)})$ opérations.

Les techniques et les résultats de Storjohann ont été étendus récemment à d'autres opérations sur les matrices polynomiales et entières, avec un accent sur le caractère déterministe des algorithmes. Ainsi, Gupta, Sarkar, Storjohann et Valeriote [Gup+12] ont donné des versions déterministes du théorème 11.4, et du calcul rapide de certaines formes échelon. Pauderis et Storjohann [PS12] ont proposé un algorithme déterministe efficace pour déterminer si une matrice entière est unimodulaire.

Néanmoins, il reste encore des problèmes ouverts. Par exemple, il n'est toujours pas connu si l'on peut calculer le polynôme minimal ou le polynôme caractéristique d'une matrice polynomiale de taille $n$ et degré au plus $d$ en $\tilde{O}(n^\theta d)$ opérations arithmétiques, que ce soit en version probabiliste ou déterministe.
Bibliographie

Dix82 DIXON, John D. (1982). « Exact solution of linear equations using *p*-adic expansions ». In: *Numerische Mathematik*, vol. 40, n°1, p. 137-141.

GJV03 GIORGI, Pascal, Claude-Pierre JEANNEROD et Gilles VILLARD (2003). « On the complexity of polynomial matrix computations ». In: *ISSAC'03: International Symposium on Symbolic and Algebraic Computation*. Éd. par J. R. SENDRA. New York: ACM Press, p. 135-142.

Gup+12 GUPTA, Somit, Soumojit SARKAR, Arne STORJOHANN et Johnny VALERIOTE (2012). « Triangular *x*-basis decompositions and derandomization of linear algebra algorithms over $K[x]$ ». In: *Journal of Symbolic Computation*, vol. 47, n°4, p. 422-453.

JV05 JEANNEROD, Claude-Pierre et Gilles VILLARD (2005). « Essentially optimal computation of the inverse of generic polynomial matrices ». In: *Journal of Complexity*, vol. 21, n°1, p. 72-86.

KV04 KALTOFEN, Erich et Gilles VILLARD (2004). « On the complexity of computing determinants ». In: *Computational Complexity*, vol. 13, n°3-4, p. 91-130.

MC79 MOENCK, Robert T. et John H. CARTER (1979). « Approximate algorithms to derive exact solutions to systems of linear equations ». In: *Symbolic and Algebraic Computation. Eurosam '79, An International Symposium on Symbolic and Algebraic Manipulation, Marseille, France, June 1979*. Vol. 72. Lecture Notes in Computer Science. London, UK: Springer-Verlag, p. 65-73.

PS12 PAUDERIS, Colton et Arne STORJOHANN (2012). « Deterministic unimodularity certification ». In: *ISSAC 2012—Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation*. ACM, New York, p. 281-288.

Sto02 STORJOHANN, Arne (2002). « High-order lifting ». In: *ISSAC'02: International Symposium on Symbolic and Algebraic Computation*. Éd. par Teo MORA. ACM Press, p. 246-254.

Sto03 — (2003). « High-order lifting and integrality certification ». In: *Journal of Symbolic Computation*, vol. 36, n°3-4, p. 613-648.

Sto05 — (2005). « The shifted number system for fast linear algebra on integer matrices ». In: *Journal of Complexity*, vol. 21, n°4, p. 609-650.

SV05 STORJOHANN, Arne et Gilles VILLARD (2005). « Computing the rank and a small nullspace basis of a polynomial matrix ». In: *ISSAC'05: International Symposium on Symbolic and Algebraic Computation*. New York: ACM Press, p. 309-316.

ZLS15 ZHOU, Wei, George LABAHN et Arne STORJOHANN (2015). « A deterministic algorithm for inverting a polynomial matrix ». In: *Journal of Complexity*, vol. 31, n°2, p. 162-173.
<figure><img src="image_32.png" /></figure>
# 12. Principe de transposition de Tellegen

## Résumé

Le principe de transposition est un ensemble de règles de transformation pour les algorithmes calculant des applications linéaires. Pour tout algorithme qui calcule des produits matrice-vecteur par une matrice fixée $\mathbf{M}$, ces règles de transformation permettent d'obtenir un algorithme dual qui calcule des produits matrice-vecteur par la matrice transposée de $\mathbf{M}$. En outre, la complexité arithmétique de l'algorithme dual est essentiellement égale à celle de l'algorithme initial.

## 12.1 Introduction

Le théorème de transposition de Tellegen affirme que, étant donnée une matrice $\mathbf{M}$ de taille $m \times n$ à coefficients dans un corps $\mathbb{K}$, sans lignes ni colonnes nulles, tout algorithme linéaire $\mathcal{A}$ qui calcule l'application linéaire $\mathbf{v} \mapsto \mathbf{M} \cdot \mathbf{v}$ de $\mathbb{K}^n \to \mathbb{K}^m$ en L opérations dans $\mathbb{K}$ peut être transformé en un algorithme dual ${}^t\mathcal{A}$ qui calcule l'application linéaire transposée $\mathbf{w} \mapsto {}^t\mathbf{M} \cdot \mathbf{w}$ de $\mathbb{K}^m \to \mathbb{K}^n$ en $L-n+m$ opérations dans $\mathbb{K}$. Ici, par algorithme linéaire on entend un algorithme qui n'utilise que des opérations linéaires en les éléments de l'entrée.

Par extension, on appelle principe de Tellegen l'ensemble des règles de transformation réalisant le passage de l'algorithme direct à l'algorithme dual.

### Motivation

Si l'algorithme utilisé pour la multiplication matrice-vecteur est l'algorithme naïf, cet énoncé devient trivial. En effet, dans ce cas :

- $\mathbf{M} \cdot \mathbf{v}$ requiert $\sum_i (\alpha_i - 1) = E - m$ opérations $\pm$ et $E$ opérations $\times$,

- ${}^t\mathbf{M} \cdot \mathbf{w}$ requiert $\sum_j (\beta_j - 1) = E - n$ opérations $\pm$ et $E$ opérations $\times$,
où $\alpha_i$ (resp. $\beta_j$) est le nombre d'éléments non nuls de la *i*-ième ligne (resp. de la *j*-ième colonne) de $\mathbf{M}$, et $\mathbf{E}$ est le nombre d'éléments non nuls de $\mathbf{M}$. Par conséquent, pour une matrice complètement générique, le théorème de Tellegen ne présente aucun intérêt. Par contre, si la matrice $\mathbf{M}$ admet une structure, il est concevable que l'on dispose d'un algorithme plus rapide que l'algorithme naïf (quadratique) pour la multiplier par un vecteur de $\mathbb{K}^n$. En utilisant le principe de transposition, on obtient aussitôt un algorithme rapide pour multiplier $^t\mathbf{M}$ par un vecteur de $\mathbb{K}^m$.

**Exemple 12.1** Considérons l'exemple suivant, qui illustre ce principe en utilisant la représentation par graphes des algorithmes linéaires. L'algorithme direct prend les valeurs $x_1$ et $x_2$ en entrée et renvoie $y_1 = ax_1 + bx_2$ et $y_2 = cx_1 + dx_2$ en sortie. Les arêtes représentent des multiplications par des valeurs scalaires $a, b, c, d$.

<figure><img src="image_33.png" /></figure>

Transposer cet algorithme revient à inverser le flot du calcul, c'est-à-dire, inverser le sens des flèches, permuter les '+' avec les '-' et les entrées avec les sorties. L'algorithme ainsi obtenu prend $y_1, y_2$ en entrée et renvoie $x_1 = ay_1 + cy_2$ et $x_2 = by_1 + dy_2$. Il calcule donc bien l'application transposée de l'application initiale. De plus, le nombre d'opérations arithmétiques utilisées par les deux algorithmes est le même : 4 multiplications et 2 additions.

<figure><img src="image_34.png" /></figure>

## Utilité

Comme nous l'avons vu au Chapitre 10, beaucoup de problèmes en calcul formel s'expriment en termes d'algèbre linéaire structurée (multiplication par un vecteur ou résolution de système). Par exemple, les opérations élémentaires sur les polynômes (multiplication, division, évaluation-interpolation, interpolation rationnelle, extrapolation, etc.) sont *linéaires* si l'on fixe l'un des opérandes : ils se codent en termes de produits matrice-vecteur $\mathbf{M} \cdot \mathbf{v}$ ou $\mathbf{M}^{-1} \cdot \mathbf{v}$, où $\mathbf{M}$ est une matrice structurée (de type Toeplitz, Hankel, compagnon, Vandermonde, Cauchy, ...)

Grâce au principe de Tellegen, comprendre, analyser et améliorer un algorithme se ramène à comprendre, analyser et améliorer son transposé. Deux sont ses utilités principales : *trouver* des solutions algorithmiques de meilleure complexité, et *clarifier* le statut de certains algorithmes existant dans la littérature, qui parfois sont simplement des transposés d'algorithmes bien connus. Cela permet ainsi le traitement algorithmique unifié des problèmes duaux. On peut résumer en disant que le principe de transposition permet de diviser par deux le nombre d'algorithmes linéaires qu'il reste à découvrir.
<figure><img src="image_35.png" /><figcaption>FIGURE 12.1 – Schéma de Horner et transposition de Tellegen : le schéma de Horner (à droite) est obtenu en transposant l’algorithme qui résout le problème dual (à gauche).</figcaption></figure>

Un exemple moins trivial

Considérons le problème de l'évaluation d'un polynôme $P \in \mathbb{K}[X]$ de degré $n$ en une valeur $a \in \mathbb{K}$. C'est une opération linéaire en les coefficients de P, de matrice $M = [1, a, ..., a^n]$ dans les bases canoniques. Le problème transposé est donc le suivant : pour une valeur donnée $x_0 \in \mathbb{K}$, calculer les produits $a^i x_0$, pour $0 \le i \le n$. Pour ce problème, un algorithme naturel consiste à multiplier $x_0$ par $a$, ensuite multiplier le résultat par $a$, et ainsi de suite.

Le transposé de cet algorithme s'obtient en parcourant l'algorithme direct en sens inverse, tout en permutant les entrées avec les sorties, et en remplaçant chaque instruction par sa *transposée*, obtenue en appliquant un nombre restreint de règles syntaxiques. Dans ce processus, les boucles *for* montantes deviennent des boucles *for* descendantes.

De cette manière, on obtient *automatiquement* un algorithme pour le problème de départ, à savoir, l'évaluation de P sur $a$. Dans notre cas, il s'avère que l'algorithme transposé coïncide avec la fameuse *méthode de Horner*, voir la Figure 12.1. Observons que l'algorithme transposé utilise $n$ opérations de plus que l'algorithme direct. Cette perte s'explique par le théorème de Tellegen : il s'agit tout simplement de la différence entre le nombre de colonnes et le nombre de lignes de $M$. Cette observation peut être utilisée pour expliquer l'optimalité de la règle de Horner.

12.2 La version en termes de graphes du principe de Tellegen

Dans cette section nous donnons une version du théorème de Tellegen dans un modèle particulier, celui des graphes de calcul (DAG).

**Définition 12.1** Un *graphe orienté acyclique*, ou DAG (de l'anglais directed acyclic graph) est un graphe orienté $G = (V, E)$ qui ne possède pas de cycle. Un DAG $\mathbb{K}$-linéaire est un DAG muni d'une fonction de poids $\lambda : E \to \mathbb{K}$.

Soit $I = \{x_1, ..., x_n\}$ l'ensemble des nœuds d'entrée de $G$. À tout sommet $v \in V$ on
associe une forme linéaire dans $K[X_1, \dots, X_n]$ de la façon suivante :

— si $v = x_i \in I$, alors $h_v := X_i$,

— si $v \in V \setminus I$, alors $h_v := \sum_{e=(w,v) \in E} \lambda(e) h_w$.

**Définition 12.2** On dit que $G$ calcule la matrice $\mathbf{M} = [a_{i,j}]$ de taille $m \times n$ si les formes linéaires associées aux nœuds de sortie sont $F_i = \sum_{j=1}^n a_{i,j} X_j$, $i = 1, \dots, m$. Le coût de $G$ est le nombre d'opérations linéaires $c(G)$ induites par $G$.

Avec ces notations, le principe de Tellegen s'énonce comme suit.

**Théorème 12.1** Soit $G$ un $K$-DAG qui calcule une matrice $\mathbf{M}$ de $M_{m,n}(K)$. Alors le graphe transposé $^tG$, obtenu en inversant le sens des flèches sans changer leur poids, calcule la matrice $^t\mathbf{M}$. De plus,

$$c(^tG) = c(G) - n + m.$$

*Démonstration.* La forme linéaire calculée par un sommet $v$ de $G$ est

$$h_v = \sum_{j=1}^{n} \left( \sum_{p \in \text{Chemin}(x_j, v)} \lambda(p) \right) X_j, \quad \text{où} \quad \lambda(p) = \prod_{e \text{ arête de } p} \lambda(e).$$

On a donc l'égalité $a_{i,j} = \sum_{p \in \text{Chemin}(x_j, F_i)} \lambda(p)$, qui montre que $^tG$ calcule bien $^t\mathbf{M}$.

L'égalité entre les coûts se déduit du fait que la quantité $c(G)-\text{card}(I(G))$ ne dépend pas de l'orientation de $G$, comme montré par le lemme suivant. ■

**Lemme 12.2** — Formule pour le coût d'un graphe de calcul. Avec les notations précédentes, on a l'égalité

$$c(G) = \text{card}(|\{e \in E | \lambda(e) \neq \pm 1\}|) + \text{card}(E) - \text{card}(V) + \text{card}(I).$$

*Démonstration.* Tout sommet $v$ qui reçoit $d(v) > 0$ arêtes contribue au coût de $G$ par $\text{card}(|\{e = (w, v) \in E | \lambda(e) \neq \pm 1\}|) + (d(v) - 1)$ opérations dans $K$.

Graphiquement, cela se voit sur la représentation de la Figure 12.2. La conclusion se déduit alors aisément de la suite d'égalités

$$\sum_{v \in V \setminus I} (d(v) - 1) = \sum_{v \in V \setminus I} d(v) - \sum_{v \in V \setminus I} 1 = \text{card}(E) - (\text{card}(V) - \text{card}(I)). \quad \blacksquare$$

## 12.3 Principe de Tellegen pour les programmes linéaires

Le but de cette section est de décrire le principe de Tellegen dans le modèle des programmes linéaires sur des machines à allocation de registres (machines à accès direct, MAD, ou random access memory, RAM, en anglais).
<figure><img src="image_36.png" /><figcaption>FIGURE 12.2 – Contribution d’une arête au coût d’un graphe de calcul.</figcaption></figure>

## Machines à registres

Commençons par définir les machines MAD à instructions linéaires.

**Définition 12.3** On fixe un entier $M$; il représente la mémoire disponible, de sorte que l'on peut accéder à $M$ registres $R_1, \dots, R_M$. On y stocke des nombres, c'est-à-dire des éléments du corps $\mathbb{K}$.

Un programme linéaire est la donnée des objets suivants :

- un sous-ensemble $R_{i_1}, \dots, R_{i_n}$ des registres que l'on appelle entrée;

- un sous-ensemble $R_{o_1}, \dots, R_{o_m}$ des registres que l'on appelle sortie;

- une suite d'instructions portant sur les registres, choisies parmi :

- $R_i = \pm R_j \pm R_k$, avec $1 \le i, j, k \le M$,

- $R_i = \lambda R_j$, avec $1 \le i, j \le M$ et $\lambda$ dans $\mathbb{K}$.

Le fonctionnement d'un tel programme est le suivant : à l'initialisation, les registres d'entrée reçoivent des valeurs $x_1, \dots, x_n \in \mathbb{K}$ et les autres sont mis à zéro. On effectue ensuite toutes les instructions, puis le programme renvoie les valeurs des registres de sortie. Noter que l'on calcule bel et bien une fonction linéaire des $x_i$.

## Transposition de programme

### Cas d'un jeu d'instructions réduit

Pour commencer, on traite le cas d'une machine avec un jeu d'instructions limité par rapport au cas général. On ne s'autorise que les instructions du type :

- $R_i = R_i \pm R_j$, avec $1 \le i, j \le M$;

- $R_i = \lambda R_i$, avec $1 \le i \le M, \lambda \in \mathbb{K}$.

Pour « transposer » un programme comportant des instructions de ce type, on interprète chacune de ces instructions comme une application linéaire $\mathbb{K}^M \to \mathbb{K}^M$. Pour motiver ce procédé, considérons un exemple. Avec $M = 3$, l'instruction $R_1 = R_1 + R_3$ s'interprète comme l'application linéaire dont la matrice est

$$ \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}; $$
$R_2$ et $R_3$ n'ont pas changé et $R_1$ devient $R_1 + R_3$. La transposée de cette application linéaire a pour matrice

$$ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix}; $$

on en déduit que l'on peut la traduire par l'instruction $R_3 = R_3 + R_1$.

De manière générale, la transposée de l'instruction $R_i = R_i + R_j$ est l'instruction $R_j = R_j + R_i$. Par écriture matricielle, on voit que l'instruction $R_i = \lambda R_i$ est inchangée par transposition, et que $R_i = R_i - R_j$ se transpose en $R_j = R_j - R_i$. On peut alors définir le programme transposé d'un programme.

**Définition 12.4** Le programme transposé $^tP$ du programme linéaire $P$ est donné par :

- les entrées de $^tP$ sont les sorties de $P$;

- les sorties de $^tP$ sont les entrées de $P$;

- les instructions de $^tP$ sont les transposées des instructions de $P$, prises en sens inverse.

La dernière condition, le retournement de la liste des instructions, traduit l'égalité matricielle $^t(AB) = ^tB \cdot ^tA$. La définition montre que si $P$ calcule une application linéaire $K^n \to K^m$, alors $^tP$ calcule bien l'application transposée $K^m \to K^n$.

**Exemple 12.2** Considérons le programme

$R_4 = R_4 + R_2$

$R_3 = 3 \times R_3$

$R_4 = R_4 + R_3$

$R_2 = 2 \times R_2$

$R_3 = R_3 + R_2$

$R_3 = R_3 + R_1$

dont les entrées sont $R_1, R_2, R_3$ et les sorties $R_3, R_4$. Ce programme calcule l'application linéaire dont la matrice est

$$ \begin{pmatrix} 1 & 2 & 3 \\ 0 & 1 & 3 \end{pmatrix}. $$

Le programme transposé s'écrit :

$R_1 = R_1 + R_3$

$R_2 = R_2 + R_3$

$R_2 = 2 \times R_2$

$R_3 = R_3 + R_4$

$R_3 = 3 \times R_3$

$R_2 = R_2 + R_4$
et on vérifie qu'il calcule l'application linéaire dont la matrice est

$$ \begin{pmatrix} 1 & 0 \\ 2 & 1 \\ 3 & 3 \end{pmatrix}. $$

### Cas général

Le cas du jeu d'instructions général se ramène au cas du jeu d'instructions réduit de manière immédiate. Ainsi, l'instruction $R_1 = R_2 + R_3$ est équivalente à la suite d'instructions

$$ \begin{align*} R_1 &= 0 \\ R_1 &= R_1 + R_2 \\ R_1 &= R_1 + R_3 \end{align*} $$

Il est donc possible de le transposer sous la forme

$$ \begin{align*} R_3 &= R_3 + R_1 \\ R_2 &= R_2 + R_1 \\ R_1 &= 0 \end{align*} $$

De même, on récrit aisément l'instruction $R_i = \lambda R_j$ en utilisant le jeu d'instructions réduit, ce qui permet de la transposer, ...

### Optimisations

Notre principe de transposition n'est pas complètement satisfaisant : au vu des règles de récriture ci-dessus, il semble que l'on puisse perdre jusqu'à un facteur 3 (en termes de nombre de lignes). On peut optimiser le code produit, afin de regagner, autant que possible, les lignes perdues. On utilise pour ce faire les règles suivantes :

**Suppression des zéros** La transformation pour passer du cas général au jeu d'instructions réduit introduit des lignes du type $R_i = 0$. On peut supprimer une telle ligne, à condition de récrire les instructions suivantes en conséquence.

**Suppression des recopies** Après avoir supprimé les zéros, il se peut qu'il reste des lignes du type $R_i = \pm R_j$. On peut également supprimer ce type de ligne, à condition de récrire de manière judicieuse les instructions qui suivent.

Muni de ces règles supplémentaires, on peut montrer qu'il est possible d'obtenir un code transposé qui fait exactement le même nombre de lignes que l'original mais cela dépasse le cadre de cet ouvrage.

## 12.4 Applications

Dans cette section, nous illustrons l'utilité pratique du principe de Tellegen à travers quelques exemples. Les techniques de transposition permettront d'engendrer des programmes linéaires pour effectuer les opérations suivantes sur les polynômes à une variable : multiplication, division euclidienne, évaluation et interpolation multipoint. La Figure 12.3 contient une liste (non exhaustive) d'opérations linéaires de base sur les polynômes à une variable et leurs problèmes transposés. Tout algorithme résolvant un problème de la colonne de gauche peut être transposé en un algorithme de même complexité pour le problème correspondant de la colonne de droite.
<table><thead><tr><th>Problème direct</th><th>Problème transposé</th></tr></thead><tbody><tr><td>multiplication mul(<span style="color:green">●</span>, <span style="color:red">●</span>)<br><span style="color:green">●</span><sup>X0</sup> × <span style="color:red">●</span><sup>Xn</sup> = <span style="color:green">●</span><sup>X0</sup> <span style="color:green">●</span><sup>Xn</sup> <span style="color:green">●</span><sup>X2n</sup></td><td>produit médian mul<sup>t</sup>(<span style="color:green">●</span>, <span style="color:red">●</span>)<br><span style="color:green">●</span><sup>X0</sup> × <span style="color:red">●</span><sup>Xn</sup> <span style="color:green">●</span><sup>X2n</sup> = <span style="color:green">●</span><sup>X0</sup> <span style="color:green">●</span><sup>Xn</sup> <span style="color:green">●</span><sup>X2n</sup> <span style="color:green">●</span><sup>X3n</sup></td></tr><tr><td>division euclidienne<br>A ↦→ A mod P</td><td>extension de récurrences<br>(a<sub>0</sub>,...,a<sub>n-1</sub>) ↦→ (a<sub>0</sub>,...,a<sub>2n-1</sub>)</td></tr><tr><td>évaluation multipoint<br>P ↦→ (P(a<sub>0</sub>),...,P(a<sub>n-1</sub>))</td><td>sommes de Newton pondérées<br>(p<sub>0</sub>,...,p<sub>n-1</sub>) ↦→ (∑ p<sub>i</sub>,...,∑ p<sub>i</sub> a<sub>i</sub><sup>n-1</sup>)</td></tr><tr><td>interpolation<br>(systèmes de Vandermonde)</td><td>décomposition en éléments simples<br>(systèmes de Vandermonde transposés)</td></tr><tr><td>décalage de polynômes<br>P(X) ↦→ P(X + 1)</td><td>évaluation de factorielles descendantes<br>P = ∑ a<sub>i</sub> X<sup>i</sup> ↦→ (P(0),...,P(n-1))</td></tr><tr><td>extrapolation sur 0,1,2,...</td><td>division modulo (X - 1)<sup>n</sup></td></tr><tr><td>q-décalage</td><td>évaluation de q-factorielles descendantes</td></tr><tr><td>composition modulaire</td><td>projection des puissances</td></tr></tbody></table>

FIGURE 12.3 – « Dictionnaire de Tellegen » pour les polynômes à une variable.

## Produit médian des polynômes

Commençons par transposer la multiplication des polynômes. Cette opération devient linéaire lorsqu'on fixe l'un des deux opérandes. Fixons une fois pour toutes un polynôme $f$ dans $K[X]$ de degré $m$. Pour $n \ge 0$ quelconque, considérons l'application de multiplication $m_{f,n}: K[X]_{\le n} \to K[X]_{\le m+n}$ qui à un polynôme $g$ associe le produit $fg$. Il est aisé de voir que la transposée de cette opération $^t m_{f,n}: K[X]_{\le m+n} \to K[X]_{\le n}$ consiste à extraire les coefficients de $X^m, X^{m+1},...,X^{m+n}$ du produit d'un polynôme de degré au plus $m+n$ par le polynôme réciproque $\tilde{f}$ de $f$.

**Exemple 12.3** Par exemple, soit $f = 2 + X + 3X^2$; la matrice de $m_{f,2}$ dans les bases canoniques de $K[X]_{\le 2}$ et $K[X]_{\le 4}$ est donnée par la matrice de type Toeplitz

$$ \begin{pmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ 3 & 1 & 2 \\ 0 & 3 & 1 \\ 0 & 0 & 3 \end{pmatrix}. $$

Sa transposée est

$$ \begin{pmatrix} 2 & 1 & 3 & 0 & 0 \\ 0 & 2 & 1 & 3 & 0 \\ 0 & 0 & 2 & 1 & 3 \end{pmatrix}. $$
Cette dernière matrice est la partie médiane de la matrice de Toeplitz suivante

$$ \begin{pmatrix} 3 & 0 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 & 0 \\ 2 & 1 & 3 & 0 & 0 \\ 0 & 2 & 1 & 3 & 0 \\ 0 & 0 & 2 & 1 & 3 \\ 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 2 \end{pmatrix}, $$

qui représente le produit du polynôme $\tilde{f} = 3 + X + 2X^2$ par un polynôme $h$ de degré 4. Cette partie médiane donne les coefficients de degré 2, 3, 4 du produit $\tilde{f}h$. Cela explique pourquoi le produit transposé est appelé *produit médian*.

Le produit médian est une opération importante, qui intervient comme brique de base dans nombre d'algorithmes sur les polynômes et les séries. Son calcul efficace est un problème important à résoudre avant de passer à la transposition d'autres algorithmes plus sophistiqués.

Par exemple, l'utilisation du produit médian permet d'améliorer, d'un facteur constant, l'efficacité de l'opérateur de Newton pour l'inverse de séries formelles, et mène aux optimisations mentionnées en page 63. Cette remarque s'applique également aux autres opérations rapides sur les séries décrites au Chapitre 3. L'accélération se répercute aussi sur la division euclidienne rapide, qui repose sur une inversion de série (Théorème 4.1). En vertu du principe de Tellegen (Théorème 12.1) appliqué aux matrices de Toeplitz, le coût du produit médian est celui du produit fg, à m opérations près. En particulier, si $m = n$, cela permet d'effectuer le calcul de la partie médiane pour le coût d'une multiplication en degré $n$, au lieu des deux qu'on attendrait naïvement. Il y a plusieurs algorithmes pour multiplier des polynômes, à chacun correspond donc un algorithme de même complexité pour calculer le produit transposé; cet algorithme transposé peut être obtenu en appliquant le principe de Tellegen, soit dans sa version graphes, soit par transformation de programmes linéaires. Plus exactement, on a le résultat général suivant.

**Théorème 12.3** Tout algorithme de complexité arithmétique $M(n)$ pour la multiplication polynomiale en degré $n$ peut être transformé en un algorithme de complexité arithmétique $M(n) + O(n)$ pour le produit médian en degrés $(n, 2n)$.

**Exercice 12.1** Vérifier le théorème dans le cas de la multiplication naïve des polynômes. Expliciter l'algorithme transposé dans ce cas. ■

**Exemple 12.4** En Figure 12.4, à gauche, est représenté un DAG qui calcule le produit de deux polynômes de degré 1 à la *Karatsuba* (en utilisant 3 multiplications au lieu de 4). Par transposition, on déduit un algorithme à la *Karatsuba* pour le produit transposé de deux polynômes de degré 1 et 2. L'algorithme direct utilise 6
<figure><img src="image_37.png" /><figcaption>FIGURE 12.4 – Produit direct et transposé à la Karatsuba, version DAG.</figcaption></figure>

opérations arithmétiques (3 additions et 3 multiplications); l'algorithme transposé en utilise 7 (4 additions et 3 multiplications). La différence de 1 opération est bien sûr prédite par le Théorème 12.1. L'algorithme de Karatsuba et son code transposé sont donnés en Figure 12.5. Tous les deux ont la même complexité $O(n^{\log 3})$.

**Exemple 12.5** La matrice de la DFT étant symétrique, à tout algorithme (linéaire) qui la calcule est associé un algorithme transposé de même coût qui la calcule également. Par exemple, à la variante de la DFT décrite au Chapitre 2 (due à Gentleman et Sande) correspond un algorithme de Cooley et Tukey. En Figure 12.6 sont représentés les DAG de ces algorithmes pour une DFT sur 4 points.

## Division avec reste et extension des récurrences à coefficients constants

Le dual du problème de la division avec reste d'un polynôme de degré N par un polynôme fixé de degré n est l'extension des récurrences linéaires à coefficients constants. Étant donnés les n premiers termes d'une suite qui vérifie une récurrence linéaire à coefficients constants d'ordre n, il s'agit de calculer les N termes suivants.

Pour une récurrence d'ordre $n = 1$, la preuve est immédiate : étendre une récurrence $u_{n+1} = au_n$, de condition initiale $u_0 = x_0$ revient à calculer la suite des valeurs $x_0, ax_0, ..., a^N x_0$ à partir de $x_0$. Or, on a vu en page 219 que le dual de ce problème est l'évaluation polynomiale en $a$, autrement dit, la division avec reste modulo $X-a$.

**Exercice 12.2** Finir la preuve dans le cas général (n quelconque). ■

**Exercice 12.3** Transposer l'algorithme rapide de division euclidienne donné dans la preuve du Théorème 4.1 (page 83). En déduire un algorithme de complexité $O(M(n))$ permettant de calculer les $2n$ premiers termes d'une suite donnée par une récurrence linéaire d'ordre $n$ à coefficients constants et $n$ conditions initiales. ■
<figure><img src="image_38.png" /><figcaption>FIGURE 12.5 – Algorithme de Karatsuba et son transposé.</figcaption></figure>

## Évaluation multipoint et interpolation

Comme vu au Chapitre 5, l'évaluation multipoint se traduit en termes matriciels par un produit entre la matrice de Vandermonde associée aux points $a_i$ et le vecteur des coefficients du polynôme P, voir la Figure 12.7. Ainsi, le problème transposé est la multiplication d'une matrice de Vandermonde transposée par un vecteur, c'est-à-dire, le calcul de sommes de Newton pondérées (on les appelle ainsi car si tous les $p_i$ valent 1, on retrouve les sommes des puissances des $a_i$). Il est possible d'exploiter algorithmiquement cette dualité, à l'aide du principe de transposition. L'idée est de proposer d'abord un algorithme rapide pour le calcul de ces sommes de Newton, et ensuite de le transposer.

Le point de départ est la série génératrice des sommes de Newton pondérées

$$ \sum_{s \ge 0} \left( \sum_{i=0}^{n-1} p_i a_i^s \right) X^{s-1} $$

dont on observe qu'elle est rationelle et vaut $Q = B/A$, où

$$ A = (X - a_0) \cdots (X - a_{n-1}) \quad \text{et} \quad B = \sum_{i=0}^{n-1} p_i \frac{A}{X - a_i} $$

D'où l'algorithme suivant : on calcule A et B et on retourne les $n$ premiers coefficients du développement en série de Taylor à l'infini de $Q = B/A$. Le calcul de A se fait en utilisant l'Algorithme 5.3 (page 100), de complexité $\frac{1}{2} M(n) \log n + O(M(n))$. Celui de B peut se faire par l'algorithme de type « diviser pour régner » représenté en en haut de la Figure 12.8, de complexité $M(n) \log n + O(M(n))$. Une solution équivalente serait de calculer simultanément A et B par l'Algorithme 5.6 en page 103. Enfin, le
<figure><img src="image_39.png" /><figcaption>FIGURE 12.6 – Dualité entre deux classes d’algorithmes pour la DFT : à gauche, l’algorithme *decimation-in-time* de Cooley–Tukey ; à droite l’algorithme *decimation-in-frequency* de Gentleman–Sande.</figcaption></figure>

$$ \begin{pmatrix} 1 & a_0 & \cdots & a_0^n \\ 1 & a_1 & \cdots & a_1^n \\ \vdots & \vdots & & \vdots \\ 1 & a_n & \cdots & a_n^n \end{pmatrix} \begin{pmatrix} p_0 \\ p_1 \\ \vdots \\ p_n \end{pmatrix} = \begin{pmatrix} P(a_0) \\ P(a_1) \\ \vdots \\ P(a_n) \end{pmatrix} \quad \text{et} \quad \begin{pmatrix} 1 & 1 & \cdots & 1 \\ a_0 & a_1 & \cdots & a_n \\ \vdots & \vdots & & \vdots \\ a_0^n & a_1^n & \cdots & a_n^n \end{pmatrix} \begin{pmatrix} p_0 \\ p_1 \\ \vdots \\ p_n \end{pmatrix} = \begin{pmatrix} \sum p_i \\ \sum a_i p_i \\ \vdots \\ \sum a_i^n p_i \end{pmatrix} $$

<caption>FIGURE 12.7 – Évaluation multipoint et sa transposée. Sommes de Newton pondérées.</caption>

développement en série de Q s'effectue en $O(M(n))$ opérations arithmétiques grâce à des itérations de Newton, comme décrit au Chapitre 3.

En résumé, on obtient un algorithme pour le calcul des sommes de Newton pondérées, de complexité

$$ \frac{3}{2} M(n) \log n + O(M(n)). $$

Par transposition, on construit un algorithme *de même complexité* pour l'évaluation multipoint. Cet algorithme est représenté graphiquement en bas de la Figure 12.8. Tout comme l'Algorithme 5.5 (page 101), le nouvel algorithme repose sur une stratégie « diviser pour régner », mais remplace, à chaque niveau de récursion, les divisions avec reste par des produits médians, moins coûteux, d'où le gain d'un facteur constant dans la complexité.

**Exercice 12.4** Écrire les détails des deux algorithmes en Figure 12.8.

<figure><img src="image_40.png" /></figure>

**Exercice 12.5** Soit $\mathbb{K}$ un corps et soit $n$ un entier. On considère le problème suivant : Étant données les valeurs en $0, 1, \dots, n-1$ d'un polynôme (inconnu) de $\mathbb{K}[X]$ de degré au plus $n$, calculer les valeurs prises par ce polynôme en $n, n+1, \dots, 2n-1$.
<figure><img src="image_41.png" /><figcaption>FIGURE 12.8 – Transposition du calcul des sommes de Newton pondérées en un calcul rapide d'évaluation multipoint : l'algorithme rapide d'évaluation multipoint par produits médians répétés (en bas) est obtenu par transposition d'un algorithme rapide pour les sommes de Newton pondérées (en haut).</figcaption></figure>

Déterminer le problème dual, donner un algorithme de complexité $O(M(n))$ pour ce dernier, et en déduire un algorithme explicite de complexité $O(M(n))$ pour le problème de départ. ■

## L'évaluation et l'interpolation ont des coûts équivalents

Un dernier exemple d'application du théorème de Tellegen, de nature plus théorique, est la preuve du fait que les problèmes d'évaluation multipoint et d'interpolation polynomiale sont équivalents du point de vue de la complexité.
**Théorème 12.4** Tout algorithme qui effectue l'évaluation multipoint (resp. l'interpolation) sur une famille fixée de points peut être transformé en un algorithme d'interpolation (resp. d'évaluation multipoint) sur la même famille de points, de même complexité à un nombre constant de multiplications polynomiales près.

Autrement dit, si $E(n)$ et $l(n)$ représentent les complexités de l'évaluation multipoint et de l'interpolation en $n$ points, alors :

$$l(n) \in O(E(n) + M(n)) \quad \text{et} \quad E(n) \in O(l(n) + M(n)).$$

*Démonstration.* Supposons qu'on dispose d'un algorithme $\mathcal{I}$ de complexité $l(n)$ pour l'interpolation sur les points $\mathbf{a} = (a_0, \dots, a_{n-1})$. Le théorème de Tellegen montre qu'on peut construire un algorithme $^{t}\mathcal{I}$ de complexité $l(n)$ pour l'interpolation transposée sur $\mathbf{a}$. On va montrer comment en déduire un algorithme $\mathcal{E}$ pour l'évaluation multipoint sur $\mathbf{a}$. L'idée est d'utiliser la factorisation matricielle

$$\mathbf{V}_{\mathbf{a}} = ({}^{t}\mathbf{V}_{\mathbf{a}})^{-1} \cdot \mathbf{H}_{\mathbf{a}}, \quad \text{où} \quad \mathbf{H}_{\mathbf{a}} = \begin{pmatrix} n & \sum_i a_i & \dots & \sum_i a_i^{n-1} \\ \sum_i a_i & \sum_i a_i^2 & \dots & \sum_i a_i^n \\ \vdots & \vdots & \ddots & \vdots \\ \sum_i a_i^{n-1} & \sum_i a_i^{n-1} & \dots & \sum_i a_i^{2(n-1)} \end{pmatrix},$$

qui suggère l'algorithmme $\mathcal{E}$ suivant :

— appliquer $\mathcal{I}$ pour calculer $A(X) = \prod_i (X - a_i)$, en complexité $l(n)$;

— calculer $\mathbf{H}_{\mathbf{a}}$ à partir des $2n-1$ premières sommes de Newton de $A(X)$, par l'algorithmme de complexité $O(M(n))$ donné en Proposition 3.10, page 67;

— calculer $\mathbf{v} = \mathbf{H}_{\mathbf{a}} \cdot \mathbf{p}$ en $O(M(n))$ opérations, comme au Chapitre 10;

— renvoyer $\mathbf{V}_{\mathbf{a}} \cdot \mathbf{p} = {}^{t}(\mathbf{V}_{\mathbf{a}}^{-1}) \cdot \mathbf{v}$ en $l(n)$ opérations en utilisant l'algorithmhe $^{t}\mathcal{I}$.

La preuve du sens inverse est laissée en exercice.

<figure><img src="image_42.png" /></figure>

**Exercice 12.6** Finir la preuve du Théorème 12.4.

<figure><img src="image_43.png" /></figure>

## Notes

La présentation de la Section 12.2 est inspirée d'un article de Kaminski, Kirkpatrick et Bshouty [KKB88]. L'algorithmhe rapide de l'Exercice 12.3 pour l'extension des récurrences est dû à Shoup [Sho91a; Sho99]. La dualité division polynomiale $\leftrightarrow$ extension de récurrences de la Section 12.4 est implicite dans l'article de Cerlienco, Mignotte, et Piras de 1987 [CMP87, §IV.3]. Elle permet de clarifier le statut de l'algorithmhe de Shoup : c'est la transposée de l'algorithmhe de Strassen [Str73] pour la division euclidienne des polynômes, fait mis en évidence par Bostan, Lecerf et Schost [BLS03].

L'algorithmhe présenté en Section 12.4, améliorant (d'un facteur constant) les algorithmes classiques d'évaluation, est dû à Bostan, Lecerf, et Schost [BLS03]. Le principe de transposition est également employé par ces auteurs pour accélérer (toujours par des facteurs constants) les algorithmes classiques pour l'interpolation et pour son problème dual (la résolution de systèmes de Vandermonde transposés). Les mêmes méthodes s'étendent au cas plus général du théorème des restes chinois [Bos+04a].
Le Théorème 12.4 est dû à Bostan et Schost [BS04]. La factorisation matricielle utilisée dans sa preuve est tirée de l'article de Canny, Kaltofen et Yagati de 1989 [CKY89].

Une autre application récente du principe de transposition est la suivante [BSS08] : si $(P_\ell)_{\ell \ge 0}$ est une famille de polynômes de $K[X]$, avec $\deg(P_\ell) = \ell$, dont la série génératrice exponentielle $\sum_{\ell \ge 0} P_\ell(X)/\ell! \cdot Y^\ell$ est de la forme<fnref n="1" /> $v(Y) \cdot \exp(X \cdot h(Y))$, alors les conversions entre la base canonique $(X^\ell)_{\ell \ge 0}$ de $K[X]$ et la base $(P_\ell)_{\ell \ge 0}$ peuvent s'effectuer en $O(M(n))$ opérations dès lors que $v$ mod $Y^n$ peut se calculer en $O(M(n))$ opérations, et pour une large classe de séries h.

## Historique

Le principe de transposition a une histoire longue et tortueuse. Cette technique de transformation des algorithmes linéaires est issue du domaine des circuits électroniques [Ant79; Bor56; PSD70] et de la théorie du contrôle [Kal59]. Le principe même remonte aux années 1950; on en trouve les traces dans un article de Tellegen [Tel52] sur les réseaux électriques. Il a été généralisé aux filtres digitaux par Fettweis [Fet71], où l'on trouve une version embryonnaire de la version par graphes. Le théorème de Tellegen a été redémontré plusieurs fois, dans divers contextes et degrés de généralité, par Fiduccia [Fid72a; Fid73], Hopcroft et Musinski [HM73], Knuth et Papadimitriou [KP81], et par Kaminski, Kirkpatrick et Bshouty [KKB88]. En calcul formel, il a été popularisé dans les travaux de Ben-Or, Tiwari [BT88], Kaltofen, Canny, Lakshman [CKY89; KCJ00], Shoup [Sho91a; Sho95; Sho99], Lecerf, Schost [LS03], Hanrot, Quercia, Zimmermann [HQZ04], Zippel [Zip90], von zur Gathen et Gerhard [GG99], ... Il est resté longtemps méconnu, comme le montre cet extrait de l'article de Wiedemann [Wie86] de 1986 : « I was not able to find an example of a linear operator that is easy to apply but whose transpose is difficult to apply ». Le nom *principe de transposition* semble avoir été introduit par Kaltofen et Shoup au début des années 1990 [Kal93b; Sho94]. Un point de vue légèrement différent sur les aspects historiques liés au principe de transposition peut être trouvé dans une note de Bernstein [Berc].

Dans la littérature du calcul formel, c'est surtout son caractère de théorème d'existence qui est exploité : connaissant un algorithme pour une certaine opération linéaire, on en déduit l'existence d'un algorithme de même complexité pour l'opération duale. Un fait mis en évidence plus récemment [BLS03] est que la transposition des programmes peut se faire de manière systématique et (quasi-)automatique. Les algorithmes sont transposés directement, d'une manière similaire à celle utilisée en différentiation automatique [GLM91], mais en tirant profit des spécificités linéaires. Par ailleurs, lorsqu'un problème linéaire doit être résolu efficacement, une démarche naturelle consiste à essayer de trouver un algorithme rapide résolvant le problème dual. Le cas échéant, une solution rapide pour le problème de départ est obtenue en re-transposant cet algorithme. C'est le point de vue adopté en Section 12.4.

## Produit médian

Le concept de *produit médian* a été introduit par Hanrot, Quercia et Zimmermann en 2004 [HQZ04] : ils ont souligné l'importance de cette *nouvelle opération* sur les polynômes. Néanmoins, en traitement du signal [Blu68; CG00], il était déjà connu

<footnote n="1">1. Une telle famille de polynômes est appelée *suite de Sheffer* [Rom84].</footnote>
que, si la multiplication est effectuée à base de DFT, le produit médian en degrés (n, 2n) a la même complexité que le produit en degré n. En effet, il se code matriciellement en un produit matrice-vecteur par une matrice de Hankel de taille n; or, celle-ci peut être plongée dans une matrice circulante de taille 2n. Cela implique que dans le modèle de multiplication polynomiale par FFT, un produit médian (n, 2n) peut être effectué en utilisant 3 DFT (deux directes et une inverse) de taille uniquement 2n (au lieu de 3n). Cette interprétation n'est plus utile dans le modèle de multiplication par l'algorithme de Karatsuba, et Hanrot, Quercia et Zimmermann [HQZ04] sont donc les premiers à avoir proposé un algorithme de type Karatsuba pour le produit médian.

## Lien avec la différentiation automatique

Il a été remarqué par Canny, Kaltofen, et Yagati [CKY89; Kal93a] que le principe de transposition est lié au mode inverse en différentiation automatique pour le calcul du gradient d'une fonction. En effet, les éléments de $^tM \cdot w$ sont les dérivées partielles de $^tV \cdot ^tM \cdot w$ par rapport aux coordonnées de $v$. Ce produit n'est autre que $^tW \cdot M \cdot v$, et le théorème de Baur-Strassen [BS83] montre que l'on peut calculer ces dérivées partielles au prix d'un surcoût linéaire.

## Lien avec la théorie de la complexité bilinéaire

Hopcroft et Musinski ont donné une version bilinéaire de l'algorithme de Tellegen [HM73] : étant donné un algorithme bilinéaire qui utilise N multiplications pour le problème, noté $(m, n, p)$, de la multiplication de deux matrices de tailles $m \times n$ et $n \times p$, il est possible d'engendrer cinq algorithmes duaux, chacun utilisant exactement N multiplications, pour les problèmes $(n, m, p)$, $(p, m, n)$, $(m, p, n)$, $(n, p, m)$ et $(p, n, m)$. Ce résultat admet la conséquence utile suivante (mentionnée en page 183) : si l'on peut calculer en N multiplications dans $\mathbb{K}$ le produit de deux matrices quelconques sur $\mathbb{K}$ de tailles $m \times n$ et $n \times p$, alors $\omega \le 3 \log_{mnp} N$.

## Dualité pour l'évaluation des monômes à plusieurs variables

Knuth et Papadimitriou [KP81] ont montré que si $M = (a_{i,j})$ est une matrice carrée inversible dont les éléments sont des entiers strictement positifs, alors, partant de $\{X_1, \dots, X_n\}$, le nombre minimum de multiplications pour évaluer les monômes

$$ \{X_1^{a_{11}} X_2^{a_{12}} \cdots X_n^{a_{1n}}, X_1^{a_{21}} X_2^{a_{22}} \cdots X_n^{a_{2n}}, \dots, X_1^{a_{n1}} X_2^{a_{n2}} \cdots X_n^{a_{nn}}\} $$

est le même que le nombre minimum de multiplications pour évaluer les monômes

$$ \{X_1^{a_{11}} X_2^{a_{21}} \cdots X_n^{a_{n1}}, X_1^{a_{12}} X_2^{a_{22}} \cdots X_n^{a_{n2}}, \dots, X_1^{a_{1n}} X_2^{a_{2n}} \cdots X_n^{a_{nn}}\}. $$

Cet énoncé peut être vu comme cas particulier du théorème de Tellegen.

## Lien entre les algorithmes de Keller-Gehrig et de Storjohann

L'itération de Keller-Gehrig (Chapitre 8, page 176) permettant le calcul de la suite

$$ s = [x, ax, a^2 x, \dots, a^N x] $$

$(N = 2^k - 1)$ est la transposée de l'algorithme « diviser pour régner » pour d'évaluation d'un polynôme $P(x) = a_0 + \cdots + a_N x^N$ en $a$ via la décomposition de type decimation-in-time $P(x) = P_0(x^2) + x \cdot P_1(x^2)$. L'algorithme « diviser pour régner » pour l'évaluation de
P en a via la décomposition de type *decimation-in-frequency* P(x) = P₀(x) + x<sup>N+1</sup>/2 P₁(x) admet comme dual l'algorithme de Storjohann (Chapitre 11, page 212) pour calculer s.

Bibliographie

<table>
  <tr>
    <td>Ant79</td>
    <td>ANTONIOU, A. (1979). Digital filters : analysis and design. McGraw-Hill Book Co.</td>
  </tr>
  <tr>
    <td>Berc</td>
    <td>BERNSTEIN, Daniel J. (2006). The transposition principle. url : http://cr.yp.to/transposition.html (visible en 2006).</td>
  </tr>
  <tr>
    <td>BLS03</td>
    <td>BOSTAN, Alin, Grégoire LECERF et Éric SCHOST (2003). « Tellegen's principle into practice ». In : ISSAC'03 : International Symposium on Symbolic and Algebraic Computation. Éd. par J. R. SENDRA. ACM Press, p. 37-44.</td>
  </tr>
  <tr>
    <td>Blu68</td>
    <td>BLUESTEIN, Leo I. (1968). « A linear filtering approach to the computation of the discrete Fourier transform ». In : IEEE Northeast Electronics Research and Engineering Meeting, vol. 10, p. 218-219.</td>
  </tr>
  <tr>
    <td>Bor56</td>
    <td>BORDEWIJK, J. L. (1956). « Inter-reciprocity applied to electrical networks ». In : Applied Scientific Research, Section B, vol. 6, p. 1-74.</td>
  </tr>
  <tr>
    <td>Bos+04a</td>
    <td>BOSTAN, A., G. LECERF, B. SALVY, É. SCHOST et B. WIEBELT (2004). « Complexity issues in bivariate polynomial factorization ». In : ISSAC'04 : International Symposium on Symbolic and Algebraic Computation. New York : ACM Press, p. 42-49.</td>
  </tr>
  <tr>
    <td>BS04</td>
    <td>BOSTAN, A. et É. SCHOST (2004). « On the complexities of multipoint evaluation and interpolation ». In : Theoretical Computer Science, vol. 329, n°1-3, p. 223-235.</td>
  </tr>
  <tr>
    <td>BS83</td>
    <td>BAUR, W. et V. STRASSEN (1983). « The complexity of partial derivatives ». In : Theoretical Computer Science, vol. 22, p. 317-330.</td>
  </tr>
  <tr>
    <td>BSS08</td>
    <td>BOSTAN, Alin, Bruno SALVY et Éric SCHOST (2008). « Power series composition and change of basis ». In : ISSAC'08 : International Symposium on Symbolic and Algebraic Computation. Éd. par David J. JEFFREY. ACM Press, p. 269-276.</td>
  </tr>
  <tr>
    <td>BT88</td>
    <td>BEN-OR, M. et P. TIWARI (1988). « A deterministic algorithm for sparse multivariate polynomial interpolation ». In : STOC'88 : ACM Symposium on Theory of Computing. ACM Press, p. 301-309.</td>
  </tr>
  <tr>
    <td>CG00</td>
    <td>CHU, Eleanor et Alan GEORGE (2000). Inside the FFT black box. Serial and parallel fast Fourier transform algorithms. CRC Press.</td>
  </tr>
  <tr>
    <td>CKY89</td>
    <td>CANNY, John F., Erich KALTOFEN et Lakshman YAGATI (1989). « Solving systems of non-linear polynomial equations faster ». In : ISSAC'89 : International Symposium on Symbolic and Algebraic Computation. ACM Press, p. 121-128.</td>
  </tr>
  <tr>
    <td>CMP87</td>
    <td>CERLIENCO, L., M. MIGNOTTE et F. PIRAS (1987). « Suites récurrentes liénaires. Propriétés algébriques et arithmétiques ». In : L'Enseignement Mathématique. II, vol. 33, p. 67-108.</td>
  </tr>
  <tr>
    <td>Fet71</td>
    <td>FETTWEIS, A. (1971). « A general theorem for signal-flow networks, with applications ». In : Archiv für Elektronik und Übertragungstechnik, vol. 25, n°12, p. 557-561.</td>
  </tr>
</table>
Fid72a FIDUCCIA, C. M. (1972). « On obtaining upper bounds on the complexity of matrix multiplication ». In : *Complexity of computer computations*. IBM Thomas J. Watson Research Center, Yorktown Heights, New York : Plenum, p. 31-40.

Fid73 — (1973). « On the algebraic complexity of matrix multiplication ». Thèse de doctorat. Brown Univ., Providence, RI, Center Comput. Inform. Sci., Div. Engin.

GG99 GATHEN, Joachim von zur et Jürgen GERHARD (1999). *Modern computer algebra*. Cambridge University Press.

GLM91 GILBERT, J.-C., G. LE Vey et J. MASSE (1991). *La différentiation automatique de fonctions représentées par des programmes*. Rapp. tech. RR INRIA 1557.

HM73 HOPCROFT, J. et J. MUSINSKI (1973). « Duality applied to the complexity of matrix multiplication and other bilinear forms ». In : *SIAM Journal on Computing*, vol. 2, p. 159-173.

HQZ04 HANROT, Guillaume, Michel QUERCIA et Paul ZIMMERMANN (2004). « The middle product algorithm I ». In : *Applicable Algebra in Engineering, Communication and Computing*, vol. 14, n°6, p. 415-438.

Kal59 KALMAN, R. E. (1959). « On the general theory of control systems ». In : *IRE Transactions on Automatic Control*, vol. 4, n°3, p. 481-491.

Kal93a KALTOFEN, Erich (1993). « Analysis of Coppersmith's block Wiedemann algorithm for the parallel solution of sparse linear systems ». In : *Applied algebra, algebraic algorithms and error-correcting codes*. Vol. 673. Lecture Notes in Computer Science. Springer-Verlag, p. 195-212.

Kal93b — (1993). « Computational differentiation and algebraic complexity theory ». In : *Workshop Report on First Theory Institute on Computational Differentiation*, p. 28-30.

KCJ00 KALTOFEN, E., R. M. CORLESS et D. J. JEFFREY (2000). « Challenges of symbolic computation : my favorite open problems ». In : *Journal of Symbolic Computation*, vol. 29, n°6, p. 891-919.

KKB88 KAMINSKI, M., D. G. KIRKPATRICK et N. H. BSHOUTY (1988). « Addition requirements for matrix and transposed matrix products ». In : *Journal of Algorithms*, vol. 9, n°3, p. 354-364.

KP81 KNUTH, Donald E. et Christos H. PAPADIMITRIOU (1981). « Duality in addition chains ». In : *Bulletin of the European Association for Theoretical Computer Science*, vol. 13, p. 2-4.

LS03 LECERF, G. et É. SCHOST (2003). « Fast multivariate power series multiplication in characteristic zero ». In : *SADIO Electronic Journal on Informatics and Operations Research*, vol. 5, n°1, p. 1-10.

PSD70 PENFIELD Jr., P., R. SPENCE et S. DUINKER (1970). *Tellegen's theorem and electrical networks*. The M.I.T. Press, Cambridge, Mass.-London.

Rom84 ROMAN, Steven (1984). *The umbral calculus*. Vol. 111. Pure and Applied Mathematics. Academic Press Inc. [Harcourt Brace Jovanovich Publishers].

Sho91a SHOUP, V. (1991). « A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic ». In : *ISSAC'91 : International Symposium on Symbolic and Algebraic Computation*. ACM Press, p. 14-21.
Sho94 SHOUP, Victor (1994). « Fast construction of irreducible polynomials over finite fields ». In : *Journal of Symbolic Computation*, vol. 17, n°5, p. 371-391.

Sho95 SHOUP, V. (1995). « A new polynomial factorization algorithm and its implementation ». In : *Journal of Symbolic Computation*, vol. 20, n°4, p. 363-397.

Sho99 — (1999). « Efficient computation of minimal polynomials in algebraic extensions of finite fields ». In : *ISSAC'99 : International Symposium on Symbolic and Algebraic Computation*. New York : ACM Press, p. 53-58.

Str73 STRASSEN, V. (1972/73). « Die Berechnungskomplexität von elementar-symmetrischen Funktionen und von Interpolationskoeffizienten ». In : *Numerische Mathematik*, vol. 20, p. 238-251.

Tel52 TELLEGEN, B. (1952). « A general network theorem, with applications ». In : *Philips Research Reports*, vol. 7, p. 259-269.

Wie86 WIEDEMANN, D. (1986). « Solving sparse linear equations over finite fields ». In : *IEEE Transactions on Information Theory*, vol. IT-32, p. 54-62.

Zip90 ZIPPEL, R. (1990). « Interpolating polynomials from their values ». In : *Journal of Symbolic Computation*, vol. 9, n°3, p. 375-403.
<figure><img src="image_44.png" /></figure>
# 13. Équations différentielles à coefficients séries

## Résumé

L'algorithmme de Newton dans un cadre différentiel permet de calculer des séries tronquées solutions d'équations différentielles en complexité quasi-optimale. Pour certaines classes particulières importantes d'équations, la complexité peut être encore abaissée.

Comme pour le Chapitre 3 sur les calculs rapides de séries, N sera utilisé pour représenter le nombre de termes à calculer, et « série » sera employé pour « série tronquée à l'ordre N ». L'efficacité des algorithmes sera mesurée par leurs complexités arithmétiques. Nous nous attacherons à expliciter la dépendance de la complexité vis-à-vis des autres paramètres (ordre *r* de l'équation, degré *d* des coefficients lorsqu'ils sont polynomiaux). Pour le cas particulier très important des équations et des systèmes différentiels linéaires, les résultats de complexité de ce chapitre sont présentés dans la Table 13.1, dans laquelle il faut comparer les complexités aux tailles des entrées et des sorties. Par comparaison, la méthode la plus directe (les coefficients indéterminés), est quadratique en N pour le cas où les coefficients sont des séries.

Tous les algorithmes sont énoncés pour des coefficients dans un anneau $\mathbb{A}$ (avec en vue le cas de coefficients polynomiaux par exemple). Dans tout ce chapitre, nous ferons l'hypothèse supplémentaire que 2, 3, ..., N sont inversibles dans $\mathbb{A}$. Pour les cas très importants d'équations différentielles linéaires à coefficients polynomiaux ou constants (et non pas séries), des algorithmes encore plus rapides sont présentés en fin de chapitre.

## 13.1 Équations différentielles d'ordre 1

Le cas de l'ordre 1 est plus simple que le cas général. Sa présentation sert d'introduction à la section suivante, où les techniques employées ici seront généralisées.
<table><thead><tr><th>Problème<br>(entrée, sortie)</th><th>coefficients<br>séries</th><th>coefficients<br>polynômes</th><th>coefficients<br>constants</th><th>taille<br>résultat</th></tr></thead><tbody><tr><td>(équation, base)</td><td>O(r<sup>2</sup>M(N))</td><td>O(dr<sup>2</sup>N)</td><td>O(rN)</td><td>rN</td></tr><tr><td>(équation, 1 solution)</td><td>O(rM(N)logN)</td><td>O(drN)</td><td>O(M(r)<sup>N</sup>/r)</td><td>N</td></tr><tr><td>(système, base)</td><td>O(r<sup>2</sup>M(N))</td><td>O(dr<sup>θ</sup>N)</td><td>O(rM(r)N)</td><td>r<sup>2</sup>N</td></tr><tr><td>(système, 1 solution)</td><td>O(r<sup>2</sup>M(N)logN)</td><td>O(dr<sup>2</sup>N)</td><td>O(M(r)N)</td><td>rN</td></tr></tbody></table>

TABLE 13.1 – Complexité de la résolution d'équations et de systèmes différentiels linéaires pour N ≫ r.

### L'équation linéaire homogène du premier ordre

Cette équation,

$$y' = A(X)y,$$

où A(X) est une série, admet la solution

$$y(X) = y(0) \exp \left( \int A(X) \right). \qquad (13.1)$$

Le calcul utilisant cette formule est dominé par celui de l'exponentielle, donc en O(M(N)) par les algorithmes du Chapitre 3.

### L'équation linéaire inhomogène du premier ordre

Il s'agit de l'équation $y' = A(X)y + B(X)$, où A et B sont des séries. La méthode de la variation de la constante consiste à poser

$$y(X) = f(X) \exp \left( \int A(X) \right),$$

et à injecter cette expression dans l'équation pour obtenir que f vérifie

$$f'(X) = B(X) \exp \left( - \int A(X) \right).$$

L'ensemble du calcul de y est donc encore borné par O(M(N)) opérations.

### Le cas non linéaire

On considère maintenant des équations non linéaires de la forme $y' = F(X, y)$, réduites dans le résultat suivant au cas linéaire ci-dessus.
**Proposition 13.1** Soit $F(X, Y) \in A[[X, Y]]$ une série à deux variables telle que, pour toute série $s(X) \in A[[X]]$, les N premiers termes de $F(X, s(X))$ et de $\frac{\partial F}{\partial y}(X, s(X))$ se calculent en $O(L(N))$ opérations. Alors, l'équation différentielle

$$y' = F(X, y), \quad y(0) = a,$$

admet une série formelle solution, et ses N premiers termes peuvent être calculés en $O(L(N) + M(N))$ opérations.

L'énoncé de cette proposition en termes d'une fonction L permet de l'adapter à différents besoins : dans le cas le pire, il est possible d'invoquer l'algorithme de Brent-Kung du Chapitre 3 pour avoir $L(N) = O(N \sqrt{N \log N} M(N))$ qui domine alors la complexité de la résolution. Cependant, pour des équations plus simples, la composition avec F et sa dérivée pourra s'effectuer en $O(M(N))$. Ce sera le cas par exemple si F s'obtient à l'aide d'un nombre constant de sommes, d'exponentielles, de logarithmes, et de puissances.

*Démonstration.* L'idée de la preuve repose sur la méthode de Newton, appliquée cette fois-ci à un opérateur :

$$\Phi : y \mapsto y' - F(X, y).$$

Le principe consiste à linéariser l'opérateur $\Phi$ au voisinage d'une approximation et à en annuler la partie linéaire, pour obtenir une nouvelle approximation. À partir de

$$\Phi(y+h) = \Phi(y) + h' - \frac{\partial F}{\partial y}(X, y)h + O(h^2),$$

il s'agit donc de calculer un incrément $h$ solution de l'équation linéaire inhomogène du premier ordre

$$h' = \frac{\partial F}{\partial y}(X, y)h - \Phi(y),$$

qui peut être résolue par la méthode de la section précédente. On déduit l'itération suivante

$$y_{k+1} := y_k + h_k, \quad h_k := \exp\left(\int \frac{\partial F}{\partial y}(X, y_k) \int (F(X, y_k) - y_k') \exp\left(-\int \frac{\partial F}{\partial y}(X, y_k)\right) \mod X^{2^{k+1}}\right).$$

Il reste à prouver la convergence quadratique. Comme pour le théorème sur l'itération de Newton sur les séries du Chapitre 3 (page 64), la preuve se fait par récurrence sur $k$ en prouvant simultanément $\Phi(y_k) = O(X^{2^k})$ et $h_k = O(X^{2^k})$. Les calculs sont les mêmes et sont donc laissés en exercice.

Pour obtenir le résultat de complexité, il suffit d'observer que l'itération requiert à chaque étape la résolution d'une équation différentielle linéaire inhomogène du premier ordre et d'invoquer le théorème « diviser pour régner ».

<figure><img src="image_45.png" /></figure>
## 13.2 Équations différentielles linéaires d'ordre supérieur et systèmes d'ordre 1

L'équation différentielle d'ordre $r$

$$a_r(X)y^{(r)}(X) + \dots + a_1(X)y'(X) + a_0(X)y(X) = 0 \quad (13.2)$$

où les coefficients $a_i$ sont des séries en X, avec $a_r(0) \neq 0$, peut être récrite comme une équation d'ordre 1

$$Y' = A(X)Y \quad (13.3)$$

en prenant pour Y le vecteur de séries $(y(X), \dots, y^{(r-1)}(X))$ et A une matrice (compagnon) de séries en X.

À l'inverse, un système (13.3) induit une relation de dépendance linéaire à coefficients des séries entre Y, Y', Y'', \dots, Y(r) si $r$ est la dimension de A (on a $r+1$ vecteurs en dimension $r$). Le vecteur Y et chacun de ses coefficients vérifient donc une équation de type (13.2).

Résoudre un problème est donc équivalent à résoudre l'autre. Dans certains cas, exploiter le caractère compagnon de la matrice induite par (13.2) mène à une meilleure complexité pour le cas des équations.

Une différence importante avec le cas des équations d'ordre 1 est que nous avons, pour les équations de type (13.2) comme pour les systèmes (13.3) deux problèmes à considérer :

- calculer une base des séries solutions;

- étant données des conditions initiales, calculer la série solution correspondante. Les complexités de ces problèmes sont liées : si l'on sait résoudre le second pour un coût C, alors on sait résoudre le premier pour un coût borné par $rC$. Si l'on sait résoudre le premier, alors une combinaison linéaire des solutions permet de résoudre le second en au plus $rN$ opérations supplémentaires.

Il est cependant utile de considérer les quatre problèmes (système ou équation, base ou solution unique) car la structure de chacun des problèmes permet de concevoir des algorithmes plus efficaces que n'en donnent les conversions directes d'un problème à l'autre.

Une méthode par « diviser pour régner »

L'équation à résoudre à précision $X^N$ est ici

$$Y' - AY = B, \quad Y(0) = v,$$

où A est une matrice de séries formelles. L'idée est de résoudre d'abord à précision moitié et de déterminer puis résoudre l'équation satisfaite par le reste. La condition initiale est d'abord traitée séparément en écrivant $Y = v + XU$ et en l'injectant dans l'équation, ce qui mène à

$$XU' + (I - XA(X))U = C, \quad C = B + A(X)v.$$

Soient maintenant $d < N$ et $U = U_0 + X^d U_1$ où $U_0$ est un polynôme de degré au plus $d-1$ en X, l'équation satisfaite par U donne :

$$XU_1' + ((d+1)I - XA(X))U_1 = -X^{-d}(XU_0' + (I - XA(X))U_0 - C).$$
**Entrée** $A_0, \dots, A_{N-p}$ dans $M_{r \times r}(\mathbb{A})$, $A = \sum A_i X^i$, $s_0, \dots, s_{N-p}$ dans $M_{r \times l}(\mathbb{A})$, $s = \sum s_i X^i$, $p \in \{1, \dots, N\}$.

**Sortie** $y$ dans $M_{r \times l}(\mathbb{A})[X]$ tel que $\deg_X y \le N - p$, et $Xy' + (pI - \text{XA})y = s + O(X^{N-p+1})$.

1. Initialiser $m$ à $N - p$.

2. Si $m = 1$ alors renvoyer $p^{-1}s(0)$ sinon :

a. Faire $d = \lfloor m/2 \rfloor$.

b. Appeler récursivement l'algorithmme avec les entrées $A$, $p+d$, $s \bmod X^d$, et mettre le résultat dans $y_0$.

c. $R = [s - Xy'_0 - (pI - \text{XA})y_0]^m_d / X^d$.

d. Appeler récursivement l'algorithmme avec les entrées $A$, $p+d$, $R$, et mettre le résultat dans $y_1$.

e. Renvoyer $y_0 + X^d y_1$.

**Algorithme 13.1** – Solution approchée de $XY' + (pI - \text{XA})Y = s$.

Si $U_0$ est une solution à précision $d$, le membre droit de l'équation est bien un polynôme. L'application de la méthodologie « diviser pour régner » amène donc naturellement à considérer l'équation plus générale

$$XY' + (pI - \text{XA})Y = R,$$

où $R$ est une série, $A$ une matrice de séries et $p \in \{1, \dots, N\}$. Celà nous conduit à l'Algorithmme 13.1, où la notation $[s]^m_d$ désigne le polynôme obtenu en ne conservant que les termes de $s$ de degrés au moins $d$ et au plus $m-1$. La suite de la résolution est détaillée dans l'Algorithmme 13.2. Les résultats de complexité se déduisent de cette méthode pour différentes valeurs de $A$ et de $\ell$. En particulier, le cas des équations $A$ une meilleure complexité que le cas général d'un système d'ordre 1 car la matrice correspondante est une matrice compagnon. Le produit d'une telle matrice par un vecteur ne coûte que $r$ opérations au lieu de $r^2$ dans le cas général.

**Théorème 13.2** — Diviser pour régner pour les équations différentielles.

Étant donnés les $N$ premiers coefficients de $A \in M_{r \times r}(\mathbb{A}[[X]])$, de $B \in M_{r \times \ell}(\mathbb{A}[[X]])$ ainsi que des conditions initiales $v \in M_{r \times \ell}(\mathbb{A})$, l'Algorithmme 13.2 calcule l'unique solution de

$$Y' - AY = B + O(X^N), \quad Y(0) = v$$

en un nombre d'opérations dans $\mathbb{A}$ borné par

1. $O(r^2 \ell M(N) \log N)$ en général;

2. $O(r \ell M(N) \log N)$ si $A$ est une matrice compagnon.

Il est possible d'améliorer les estimations de complexité pour cet algorithme pour la valeur $\ell = r$ (calcul de toute une base des solutions), mais dans ce cas, il vaut mieux employer les algorithmes de la section suivante, qui sont plus efficaces.

**Démonstration.** Il suffit de prouver le résultat pour $N$ une puissance de 2, et le cas
**Entrée** $Y_0, A_0, \dots, A_{N-2}$ dans $M_{r \times r}(\mathbb{A})$, $A = \sum A_i X^i$.

**Sortie** $Y = \sum_{i=0}^{N-1} Y_i X^i$ dans $M_{r \times r}(\mathbb{A})[X]$ tel que
$Y' = AY + O(X^{N-1})$, et $Z = Y^{-1} + O(X^{N/2})$.

$$1. Y = (I + XA_0)Y_0.$$

$$2. Z = Y_0^{-1}.$$

$$3. m = 2.$$

4. Tant que $m \le N/2$ faire :

$$Z = Z + Z(I_r - YZ) \mod X^m;$$

$$Y = Y - Y\left(\int Z(Y' - A \mod X^{2m-1}Y) \right) \mod X^{2m};$$

$$m = 2m.$$

5. Renvoyer Y, Z.

**Algorithme 13.3** – Résolution de $Y' = A(X)Y$, $Y(0) = Y_0$ par itération de Newton.

Cette méthode requiert le calcul de l'inverse d'une base de solution, inverse qui peut être calculé simultanément par l'itération de Newton :

$$Z_{k+1} = Z_k + Z_k(I - Y_k Z_k).$$

Une version précise de la méthode, avec les ordres de troncature, est donnée dans l'Algorithme 13.3.

**Lemme 13.3** — Correction de l'algorithme. Soit $m$ un entier pair, soient $y$ et $z$ dans $M_{r \times r}(\mathbb{A}[X])$ tels que

$$I - yz = O(X^{m/2}), \quad y' - Ay = O(X^m).$$

Soient ensuite $Y$ et $Z$ définis par

$$Z := z(2I - yz) \mod X^m, \quad Y := y\left(I - \int Z(y' - Ay) \right) \mod X^{2m}.$$

Alors $Y$ et $Z$ vérifient

$$I - YZ = O(X^m), \quad Y' - AY = O(X^{2m-1}).$$

*Démonstration.* D'après l'hypothèse sur $y' - Ay = O(X^m)$, $Y = y + O(X^m)$ et donc la convergence de l'inversion est donnée par

$$\begin{align*} I - YZ &= I - y(z + z(I - yz)) + O(X^m) \\ &= (I - yz) - yz(I - yz) + O(X^m) \\ &= (I - yz)^2 + O(X^m) = O(X^m). \end{align*}$$

La seconde propriété s'obtient en injectant la définition de $Y$ dans $Y' - AY$ (on pose
$Q = \int Z(y' - Ay)) :$

$$
\begin{align*}
Y' - AY &= y' + y'Q - Ay - AyQ - yZ(y' - Ay) + O(X^{2m}), \\
&= (I - yZ)(y' - Ay) - (y' - Ay)Q + O(X^{2m}), \\
&= O(X^m)O(X^m) + O(X^m)O(X^{m-1}) + O(X^{2m}) = O(X^{2m-1}). \quad \blacksquare
\end{align*}
$$

Comme d'habitude, l'application du théorème « diviser pour régner » mène au résultat
de complexité.

**Théorème 13.4 — Newton pour les systèmes différentiels linéaires.** Les N premiers termes d'une base de solutions de $Y' = AY$ sont calculés en $O(MM(r, N))$ opérations dans $\mathbb{A}$ par l'Algorithmme 13.3.

La notation MM(r, N) représente la complexité du produit de matrices r × r de poly-
nômes de degré au plus N ; des bornes non triviales sont données au Chapitre 5 :

$$
\mathrm{MM}(r, N) = O(r^{\theta}N + r^{2}M(N)).
$$

Exercice 13.1 L'exponentielle d'une série est obtenue en O(M(N)) opérations par cet algorithme lorsque r = 1. Vérifier que la constante est meilleure que celle de l'algorithme du Chapitre 3.

Exercice 13.2 Le cas d'un système inhomogène $Y' = AY + B$, où $B$ est une matrice de séries, s'obtient à partir de l'algorithme précédent par variation de la constante. Étudier les précisions requises dans les troncatures, et donner le résultat de complexité.

13.3 Cas particuliers

Équations différentielles linéaires à coefficients polynomiaux

Pour ces équations, la méthode des coefficients indéterminés donne une complexité
linéaire en N. L'explication tient à une propriété vue au Chapitre 14 : les solutions
séries de ces équations ont des coefficients qui vérifient des récurrences linéaires. Les
résultats du Chapitre 15 sur les récurrences linéaires à coefficients polynomiaux s'ap-
pliquent alors, et en particulier les N premiers coefficients d'une solution s'obtiennent
en O(drN) opérations si d est une borne sur le degré des $q_i$. Le Chapitre 16 revient
sur ces questions.

Le cas matriciel se traite de la même façon, la récurrence ci-dessus étant alors
à coefficients matriciels. Les complexités sont de même nature, avec en outre la
possibilité d'utiliser un produit rapide de matrices dans le cas du calcul d'une base
de solutions. Les résultats correspondants sont donnés dans la Table 13.1.