Cours 2: FFT

Convolution

On définit la convolution de deux séquences de valeurs:

A = (a_i)_{i ∈ [m;n]}
B = (b_j)_{j ∈ [p;q]}
A*B = (Σ_{i+j = k} a_i b_j)_{k ∈ [m+q;n+q]}

Visuellement on combine les coeffs de A et B sur les diagonales d’une matrice.

Application 1 Produit de polynomes.

Application 2 Lissage d’une série de mesures.

• A est une série de mesures, imprécises ou bruitées, qu’on souhaite lisser.

• Le filtre (ou flou) Gaussien de A est défini par:

  a'_i = 1/Z Σ_{δ ∈ [-k;k]} a_i e^{-δ²}

• En fait c’est A*B où B est le masque défini par:

  b_δ = 1/Z e^{-δ²}

Application 3 Combinaison d’histogrammes. - On prend a_i le nombre de solutions valant i points à l’exo A. - On prend b_j le nombre de solutions valant j points à l’exo B. - Le nombre c_k de solutions combinées à l’exo A et B valant k points est donné par la convolution A*B.

Le calcul naïf de la convolution est quadratique. On va voir un algo magique en O(n log n): la FFT.

Multiplication d’entiers

Multiplication de deux entiers donnés via leurs n bits. De nouveau, le calcul naïf est quadratique.

Karatsuba propose en 1960 une méthode diviser-pour-régner:

x = x_0 2^n + x_1
y = y_0 2^n + y_1
xy = x_0 y_0 2^{2n} + (x_0 y_1 + x_1 y_0) 2^n + x_1 y_1

Si on calcule le terme du milieu via deux multiplications on retrouve O(n²), mais il s’exprime aussi comme (x_0 + x_1) (y_0 + y_1) - x_0 y_0 - x_1 y_1 i.e. seulement une multiplication de plus. Si on fait l’analyse de complexité on a un T(n) = 3 T(n/2) + O(n) soit O(n ^ log_2 3) ou encore O(n ^ 1.59).

On détaille le calcul:

• au pire T(n) opérations avec T(n) = 3 T(n/2) + O(n)
• T(2^k) = Σ_{i ∈ [0;k]} 3^i O(2^k-i)
• T(2^k) = O(2^k Σ_i (3/2)^i)
• T(2^k) = O(2^k {(3/2)^{k+1} - 1}/{3/2 - 1})
• T(2^k) = O(3^k) = O((2^k)^log_2 3)

L’histoire ne s’arrête pas là:

• 1968: Schönhage & Strassen en O(n * log n * log log n)
• 2019: Harvey & van der Hoeven en O(n * log n)

Les deux algorithmes s’appuient sur la FFT. Leur intérêt est purement théorique, par contraste avec l’utilisation pratique intensive de la FFT en traitement du signal.

L’existence d’un algo en O(n * log n) avait été conjecturée en 1971 par Schönhage et Strassen. On ne sait pas si c’est optimal.

FFT

Vers la FFT

On voit la convolution comme une multiplication de polynômes. Considérons pour simplifier des polynomes sur les réels. Pour la complexité on compte en nombre d’opérations sur les réels: c’est abusé mais en fait ça ne se passe pas très différemment si on n’est pas sur les réels – et en traitement du signal, on implémente les réels par des flottants donc c’est OK de voir les opérations dessus comme primitives, et on s’assied sur les problèmes d’imprécision/arrondis.

Schéma général pour deux polynomes de taille N (degré N-1):

• choisir 2N valeurs
• évaluer chaque polynome sur ces valeurs
• multiplier terme à terme
• interpoler pour obtenir les coefficients du polynome produit

On n’a pas forcément avancé car on ne sait pas faire l’évaluation en temps mieux que quadratique… pareil pour l’interpolation. Mais on va y arriver en choisissant bien les 2N valeurs et en utilisant l’approche diviser-pour régner.

Cela va être utile de passer dans les complexes…

La FFT

Transformée de Fourier discrète (DFT)

Pour un polynome A de longueur L et ω une racine n-ème de l’unité:

DFT_ω(A) = (A(ω^0),A(ω^1)...A(ω_L))

FFT: calcul de la DFT en O(L log L)

Supposons L pair: L = 2N.

On choisit d’évaluer sur les 2N racines 2N-èmes de l’unité.

Ω_p = { ω_{i,p} | 0 ≤ i < p }
ω_{i,p} = e^{2 i j π / p}

On décompose A en ses parties paire et impaire:

A(X) = A_p(X^2) + X A_i(X^2)
A_p(X) = a_0 + a_2 X + a_4 X^2 + a_6 X^3 + ... a_{2N-2} X^N
A_i(X) = a_1 + a_3 X + a_5 X^2 + a_7 X^3 + ... a_{2N-1} X^N

On récupère les valeurs v_j et v'_j de A_p et A_i sur les racines N-èmes de l’unité. On peut trouver les valeurs des A(ω_{j,2N}) à partir de là:

A(ω_{j,2N})  = A_p(ω_{j,N}) + ω_{j,2N} A_i(ω_{j,N})

Plus concrètement, à partir des v_j,v'_j pour 0≤j<N:

A(ω_{j,2N})   = v_j + ω_{j,2N} v'_j
A(ω_{j+N,2N}) = v_j + ω_{j+N,2N} v'_j
              = v_j - ω_{j,2N} v'_j

On a exploité le fait que A_p(ω_{j+N,N}) = A_p(ω_{j,N}) = v_j et de même pour A_i, puis ω_{j+N,2N} = -ω_{j,2N}.

On peut écrire tout ça en pseudo-code simple, avec ω = ω_{1,2N}:

α := 1
pour j = 0 à N-1:
  y_j := v_j + α v'_j
  y_{j+N} := v_j - α v'_j
  α := α * ω

Interpolation

A partir des valeurs v_j de la DFT(C) on cherche à trouver la représentation de C par ses coefficients.

La remarque magique (qui demande d’avoir compris qu’on peut faire la transformée de Fourier sur les puissances de n’importe quelle racine L-ème de l’unité, et pas forcément les puissances de ω_{1,L}):

Evaluation(A) = DFT_ω(A)
Interpolation(V) = 1/L DFT_{1/ω}(V)

Posons en effet:

D(X) = Σ_{j∈[0;L-1]} v_j X^j
D(X) = Σ_{j∈[0;L-1]} C(ω_{j,L}) X^j

Vérifions:

D(ω_{k,L}) = Σ_j C(ω_{j,L}) ω_{k,L}^j
D(ω_{k,L}) = Σ_j C(ω_{j,L}) ω_{j k,L}
D(ω_{k,L}) = Σ_j (Σ_l c_l ω_{j l, L}) ω_{j k,L}
D(ω_{k,L}) = Σ_l c_l (Σ_j ω_{j l, L} ω_{j k,L})
D(ω_{k,L}) = Σ_l c_l (Σ_j ω_{j(k+l), L})

La somme sur j vaut L si ω_{k+l,L} = 1, càd si l = L-k. Pour les autres valeurs de l, comme ω_{_,L} est une racine N-ème:

Σ_j ω_{j(k+l), L} = (ω_{k+l,L}^L - 1) / (ω_{k+l,L} - 1) = 0

Donc:

D(ω_{k,L}) = c_{L-k} L
c_k = D(ω_{L-k,L}) / L

Pour les amateurs, on peut reformuler tout ça en termes de matrice de Van der Monde.

Références

Le livre Algorithm Design de Kleinberg & Tardos.

Pour aller plus loin

On peut réaliser la FFT avec un circuit: ça donne une implémentation en place et parallèle.

Implémentation sans complexes ni réels: se placer dans un corps Z/nZ avec n bien choisi. C’est la base de l’algo de Schönhage-Strasse (et, je crois, Harvey-Hoeven).

Liens avec le traitement du signal, les maths, etc. https://en.wikipedia.org/wiki/Fast_Fourier_transform