Transformée de Fourier discrète et multiplications

27 décembre 2022 — FFT, Python

Dans cet article, nous allons introduire la transformée de Fourier discrète, qui associe à une séquence de nombres complexes une autre séquence de complexes de même longueur $n$.

L’algorithme de transformation de Fourier rapide (FFT) permet de calculer cette transformée en $O(n \log n)$. En particulier, il peut être utilisé pour multiplier rapidement des polynômes, ou des entiers.

Le code source des fonctions présentées est disponible ici.

Transformée de Fourier discrète

La transformation de Fourier discrète associe à une séquence $(a_0, ..., a_{n-1})$ de nombres complexes sa transformée $(A_0, ..., A_{n-1})$.

Pour ce faire, on considère le polynôme ${P = a_0 + a_1 X + a_2 X^2 + ... + a_{n-1} X^{n-1}}$ associé à la séquence initiale et on évalue ce polynôme en $n$ points bien choisis pour obtenir la transformée.

Plus précisément, on choisit les $n$ racines $n$-èmes de l’unité ${1, \omega, \omega^2, ... \omega^{n-1}}$.

On définit ici $\omega$ comme une racine $n$-ème primitive de l’unité, par exemple ${\omega = e^{-\frac{2i\pi}{n}}}$.

On a donc : ${A_k = P(\omega^k) = \displaystyle\sum_{m=0}^{n-1} a_m e^{-2i\pi\frac{mk}{n}}}$

Cette relation s’écrit aussi matriciellement :

La transformation de Fourier inverse consiste à inverser cette matrice, je vous laisse vérifier que la matrice ci-dessous convient :

Par exemple, dans le cas $n = 2$ la transformée de Fourier discrète de $(a, b)$ est $(A, B) = (a + b, a - b)$ constituée de la somme et de la différence. La transformée inverse s’obtient par demi-somme et demi-différence : ${(a, b) = (\frac{A + B}{2}, \frac{A - B}{2})}$.

Produit de polynômes

Si $P$ et $Q$ sont deux polynômes de degrés inférieurs à $d_1$ et $d_2$, leur produit est un polynôme de degré inférieur à $d_1 + d_2$, dont au plus les $n = d_1 + d_2 + 1$ premiers coefficients sont non nuls.

La transformée de Fourier discrète (de taille $n$) de $PQ$ s’obtient en évaluant les $(PQ)(\omega^k) = P(\omega^k)Q(\omega^k)$ pour $0 \leq k \lt n$.

En particulier, on remarque que : la transformée de $PQ$ est égale à la multiplication point par point des transformées de $P$ et de $Q$.

Notons que toutes les transformées considérées ici sont de taille $n$, une séquence de taille $m < n$ doit donc être complétée par $n - m$ zéros pour calculer sa transformée.

La transformée de Fourier rapide

On suppose ici que $n$ est une puissance de $2$ : $n = 2^l$ pour un certain $l \in \N$.

On peut écrire

où ${P_0 = a_0 + a_2 X + ... + a_{n-2} X^{\frac{n}{2}-1}}$ et ${P_1 = a_1 + a_3 X + ... + a_{n-1} X^{\frac{n}{2} - 1}}$.

En particulier, on a pour $0 \leq k \lt n$ : $P(\omega^k) = {P_0(\omega^{2k}) + \omega^k P_1(\omega^{2k})}$

Dans l’expression précédente, $\omega^{2k} = (\omega^2)^k$ où $\omega^2$ est une racine $\frac{n}{2}$-ème primitive de l’unité.

On a donc : $P(\omega^k) = {P_0((\omega^2)^{k \mod \frac{n}{2}}) + \omega^k P_1((\omega^2)^{k \mod \frac{n}{2}})}$

Cette expression fait intervenir l’évaluation des polynômes $P_0$ et $P_1$ sur $\frac{n}{2}$ racines $\frac{n}{2}$-èmes de l’unité.

C’est précisément ce que calculent les transformées de taille $\frac{n}{2}$ de $P_0$ et $P_1$. On en déduit donc un algorithme de calcul de la transformée qui réduit le problème en deux sous-problèmes de taille $\frac{n}{2}$ :

from math import e, pi
def fft(L, sign=-1):
    n = len(L)
    if n == 1:
        return L
    w_n = e**(sign * 2j * pi / n)
    w = 1
    F_even = fft(L[::2], sign)
    F_odd = fft(L[1::2], sign)
    F = [0] * n
    for k in range(n // 2):
        F[k] = F_even[k] + w * F_odd[k]
        F[k + n // 2] = F_even[k] - w * F_odd[k]
        w = w * w_n
    return F

On désignera maintenant par FFT la transformée de Fourier discrète d’une séquence obtenue par cette méthode.

Note : l’algorithme implémenté ici est celui décrit dans le chapitre 30 de Introduction to Algorithms, que je vous recommande pour une présentation plus détaillée.

Complexité de l’algorithme de FFT

La fonction fft effectue un nombre proportionnel à $n$ d’opérations et exécute deux appels récursifs à elle-même sur des séquences de taille $\frac{n}{2}$.

Le nombre d’opération effectuée est donc de l’ordre de ${n + 2 \frac{n}{2} + 2^2 \frac{n}{2^2} + ... + 2^l \frac{n}{2^l}} = n(l+1) = O(n \log n)$

Transformation inverse

On a vu plus haut que la matrice d’inversion s’obtenait en considérant la matrice de la transformée, en remplaçant $\omega$ par son conjugué et en divisant par $n$.

On en déduit la fonction inverse de la transformée ci-dessus :

def ifft(L):
    return [c / len(L) for c in fft(L, sign=1)]

Produit de polynômes

En particulier, on en déduit un algorithme qui calcule le produit de deux polynômes en procédant ainsi :

1. Déterminer le plus petit $n$ supérieur à $\deg P + \deg Q + 1$ qui soit une puissance de $2$ (en $O(\log n)$ instructions).
2. Calculer la FFT de $P$ et $Q$, en ajoutant des zéros de padding pour avoir des séquences de tailles $n$ (en $O(n \log n)$ instructions).
3. Multiplier point par point les deux FFT obtenues (en $O(n)$ instructions).
4. Prendre l’inverse de cette FFT : on obtient les coefficients du polynôme $PQ$ (en $O(n \log n)$ instructions).

Cet algorithme a une complexité en $O(n \log n)$ alors que la multiplication naïve (en appliquant la formule) amènerait à une multiplication en $O(n^2)$ opérations.

En voici une implémentation en Python :

def multiply(P, Q):
    n = 1
    while n < len(P) + len(Q) - 1:
        n <<= 1
    P = P + [0] * (n - len(P))
    Q = Q + [0] * (n - len(Q))
    F = [a * b for a, b in zip(fft(P), fft(Q))]
    return ifft(F)

>>> multiply([1, 2, 3], [4, 5])
[(4-6.661338147750939e-16j),
 (13+6.123233995736765e-17j),
 (22+6.661338147750939e-16j),
 (15-6.123233995736765e-17j)]

Notons que cet algorithme a le défaut de propager des erreurs d’approximation des flottants comme ci-dessus. On peut nettoyer les coefficients obtenus à l’aide de Numpy :

>>> import numpy as np
>>> P = multiply([1, 2, 3], [4, 5])
>>> np.trim_zeros(np.round(np.real(P)).astype(np.int32), "b")
array([ 4, 13, 22, 15], dtype=int32)

Algorithme de Schönhage-Strassen

Une variante de l’algorithme décrit précédemment est utilisée dans un algorithme de multiplication d’entiers en ${O(n \log n \cdot \log \log n)}$, l’algorithme de Schönhage-Strassen.

Nous avons vu précédemment que le calcul dans $\Complex$ de la transformée introduit inévitablement des erreurs d’approximations des flottants. Ceci n’est pas acceptable pour un algorithme de multiplication d’entiers qui doit obtenir un résultat exact.

Transformée sur un anneau quelconque

La transformée que nous avons décrite repose sur des opérations d’additions et multiplications et sur l’existence de racines primitives $n$-èmes de l’unité.

Les résultats que nous avons obtenus, notamment la formule d’inversion et l’algorithme de calcul du produit de polynômes restent vrais en utilisant des séquences à valeurs dans d’autres anneaux que $\Complex$.

Il est notamment possible de considérer des séquences dans des anneaux finis, où les calculs peuvent être effectués de façon exacte par un ordinateur.

C’est ce type d’anneau qui est utilisé dans l’algorithme de Schönhage-Strassen (voir ici).

Notons notamment l’utilisation de cet algorithme dans V8 (le moteur Javascript utilisé notamment par Google Chrome et Node.js) pour multiplier des entiers de plus de $1500$ chiffres (voir le code source).

Maths et applications, avec les mains et avec du code 💻
Suivre sur Twitter

• ← La géométrie projective, d'Urbino à OpenGL
• Compression JPEG et composantes principales →