Convolution, transformée de Fourier rapide et polynômes (2022)

 (alvarorevuelta.com)
1 points par GN⁺ 2024-07-02 | 1 commentaires | Partager sur WhatsApp
  • Développer un polynôme de grand degré à la manière apprise au lycée impose de multiplier toutes les paires de termes, ce qui fait vite de O(n²) un goulot d’étranglement
  • La multiplication des vecteurs de coefficients de polynômes est équivalente à la convolution de signaux discrets, et le résultat de [2, 3, 4] et [5, 6, 7] est [10, 27, 52, 45, 28]
  • La DFT projette un signal discret dans le domaine fréquentiel, et la FFT calcule cette même transformation en O(n log n), ce qui change tout pour les grandes entrées
  • Comme la convolution dans le domaine temporel devient une multiplication terme à terme dans le domaine fréquentiel, on peut multiplier des polynômes plus vite en passant par une FFT, en multipliant, puis en revenant avec une IFFT
  • Pour les petits degrés, le coût aller-retour FFT/IFFT peut annuler le gain, mais plus le degré augmente, plus l’approche FFT devient efficace

Pourquoi la multiplication de polynômes devient lente

  • Un polynôme P(x) s’exprime comme une somme de termes formés d’un coefficient a_k et d’une puissance de la variable x
    • Exemple : P(x)=5x²+2x+9 est un polynôme de degré 2
    • Selon la convention d’écriture, le vecteur de coefficients peut s’écrire [5, 2, 9] ou [9, 2, 5]
  • L’addition et la soustraction sont relativement simples, car il suffit d’additionner ou de soustraire les termes de même degré
    • En Python, on peut parcourir les coefficients avec zip(p, q) et calculer a + b ou a - b
    • Si les degrés diffèrent, on peut utiliser zip_longest
  • La multiplication, elle, exige de multiplier chaque terme par tous les autres puis de regrouper les termes de même degré, ce qui augmente fortement le volume de calcul
    • (2x²+3x+4) × (5x²+6x+7) donne 10x⁴+27x³+52x²+45x+28
    • La complexité de cette méthode est O(n²), et le nombre de multiplications nécessaires augmente avec le degré

Vecteurs de coefficients et convolution

  • Dans le domaine discret, la convolution de deux signaux p et q est définie par y[n]=Σ p[k]·q[n-k]
  • Le calcul consiste à inverser q, puis à le faire glisser de gauche à droite sur p en additionnant les produits des éléments qui se chevauchent
  • Les signaux de l’exemple sont les suivants
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • En inversant puis en décalant q, on obtient chaque coefficient de sortie dans l’ordre suivant
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • Le résultat de la convolution est y = [10, 27, 52, 45, 28]
    • Il correspond aux coefficients de 10x⁴+27x³+52x²+45x+28, obtenu par multiplication de polynômes
    • On peut donc voir la multiplication de polynômes comme une convolution de vecteurs de coefficients

Transformée de Fourier et FFT

  • La transformée de Fourier convertit un signal du domaine temporel vers le domaine fréquentiel
    • Dans une perspective temporelle, on regarde le signal comme une suite de valeurs à différents instants
    • Dans une perspective fréquentielle, on l’interprète comme une somme d’oscillations de fréquences différentes
  • Les fréquences d’oscillation s’expriment avec des sinus et des cosinus, chacun avec son coefficient et sa phase
  • Si l’on applique une FFT à une onde sinusoïdale pure de 5 Hz, elle apparaît comme un pic de type delta à 5 Hz dans le domaine fréquentiel
    • Cela montre qu’une sinusoïde dans le domaine temporel peut être représentée par une unique composante sinusoïdale de 5 Hz
  • La terminologie associée se distingue ainsi
    • Fourier Transform (FT) : transformée de Fourier définie dans le domaine continu
    • Discrete Fourier Transform (DFT) : transformée de Fourier définie pour les signaux discrets
    • Fast Fourier Transform (FFT) : algorithme qui calcule la DFT en O(n log n) au lieu de O(n²)
  • La DFT transforme un signal discret dans le temps x[n] en sa représentation fréquentielle X[k]
    • Chaque X[k] est calculé en multipliant les échantillons d’entrée par des nombres complexes représentant une fréquence donnée, puis en les additionnant

Passer à une multiplication dans le domaine fréquentiel

  • Le principal avantage de la DFT et du domaine fréquentiel est qu’ils permettent de transformer une convolution en multiplication terme à terme
    • Convoluer deux signaux dans le domaine temporel revient à multiplier ces deux signaux dans le domaine fréquentiel
    • Une multiplication se calcule plus rapidement qu’une convolution
  • La procédure pour accélérer la multiplication de polynômes est la suivante
    • Transformer les polynômes dans le domaine fréquentiel avec une FFT : O(n log n)
    • Multiplier terme à terme dans le domaine fréquentiel : O(n)
    • Revenir dans le domaine temporel avec une IFFT : O(n log n)
  • Globalement, utiliser la FFT permet d’effectuer la multiplication de polynômes avec une complexité de O(n log n)
  • Pour les grands polynômes, c’est plus rapide que la multiplication scolaire en O(n²)

Implémentation Python et benchmark

  • multiply_naive utilise une double boucle pour multiplier toutes les paires de coefficients et ajouter le résultat à la position i + j
    • La longueur du résultat est len(p) + len(q) - 1
    • Sa complexité est O(n²)
  • multiply_fft effectue la multiplication des coefficients à l’aide de la FFT/IFFT
    • Il calcule une longueur qui est une puissance de 2 et qui vaut au moins len(p) + len(q) - 1, afin de pouvoir contenir le résultat
    • Il applique un padding aux deux entrées avec np.pad
    • Il multiplie terme à terme les valeurs transformées par np.fft.fft
    • Il revient ensuite avec np.fft.ifft, puis arrondit la partie réelle pour retrouver des coefficients entiers
  • Avec l’exemple p = [2, 3, 4], q = [5, 6, 7], les deux méthodes renvoient [10, 27, 52, 45, 28]
  • Dans le benchmark, la comparaison porte entre la méthode FFT et multiply_convolve, qui utilise np.convolve, plutôt qu’avec multiply_naive
    • multiply_naive repose sur des boucles Python, trop lentes pour être comparées directement à une approche FFT basée sur np.fft.fft
    • np.convolve effectue la même opération dans du code C bas niveau
  • Le degré varie dans la plage range(1, 30000, 1000), et pour chaque degré, deux polynômes sont générés avec des coefficients aléatoires compris entre 1 et 999999
    • Chaque méthode est mesurée sur une moyenne de n_runs = 5
    • Pour les faibles degrés, le coût du double passage FFT/IFFT peut empêcher l’approche FFT d’être avantageuse
    • Quand le degré augmente, l’approche FFT devient nettement plus efficace

À lire aussi

1 commentaires

 
GN⁺ 2024-07-02
Avis de Hacker News

  • Ce qui me gêne toujours dans ce genre d’explications, c’est qu’elles oublient généralement les erreurs numériques.
    On ne peut pas simplement abstraire la multiplication des coefficients comme étant en « temps constant ». Si l’on fait ça, on peut tout aussi bien abstraire toute la multiplication dès le départ. En tenant compte de la précision numérique, on est plutôt proche de O(n (log n)^3) [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • La borne d’erreur citée dans cet article est beaucoup trop pessimiste. La dernière édition de Knuth contient la bonne borne, parce que je la lui ai signalée
    • Ce serait bien si les opérations à base de quaternions présentées dans l’article de l’OP pouvaient permettre de réduire, voire d’éliminer, les erreurs de multiplication [1],[2],[3]
      [1] One-Dimensional Quaternion Discrete Fourier Transform and an Approach to Its Fast Computation:
      https://www.mdpi.com/2079-9292/12/24/4974
      [2] Convolution Theorems for Quaternion Fourier Transform: Properties and Applications:
      https://onlinelibrary.wiley.com/doi/10.1155/2013/162769
      [3] On the Matrix Form of the Quaternion Fourier Transform and Quaternion Convolution:
      https://arxiv.org/abs/2307.01836
    • Si les coefficients sont des entiers, on peut obtenir un résultat exact avec une NTT en utilisant un modulus suffisamment grand, et le temps de multiplication peut même être plus rapide, surtout en hardware
    • C’est pour ça qu’on distingue informatique théorique et génie logiciel :)

  • On peut multiplier de longs nombres entre eux avec cette méthode. Le point clé est que la multiplication de polynômes est identique à la multiplication longue ordinaire de grands nombres, mais sans retenues (carry).
    Par exemple, si l’on a un nombre à 1000 chiffres, on prend chaque chiffre comme coefficient d’un polynôme à 1000 éléments. Ensuite, on peut multiplier ces polynômes avec la méthode FFT décrite dans l’article. Pour reconvertir le résultat en nombre, il faut gérer les retenues. Si un élément dépasse 10, on reporte l’excédent sur le chiffre suivant, puis on transforme les coefficients en nombre.
    L’idée de base est celle-ci, mais il y a des subtilités concernant la précision nécessaire pour les retenues et la garantie que l’arrondi du résultat de la FFT à l’entier le plus proche donne bien le bon résultat. C’est ainsi que GMP, la bibliothèque de référence du domaine, effectue les multiplications de grands nombres

    • Comme tu le dis, on peut représenter un nombre décimal comme un polynôme avec x=10, donc ça se comprend. Par exemple, 983 = 9x^2 + 8x + 3, soit [9, 8, 3].
      Je me demande à partir de quelle taille de nombres cela devient réellement pertinent, et dans quels usages

  • Si vous ne l’avez pas encore vue, cette vidéo vaut le détour
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    Elle dérive l’algorithme FFT à partir de la multiplication de polynômes, et c’est vraiment excellent. Je la revois à peu près tous les six mois

  • La propriété de la FFT selon laquelle « la convolution est une multiplication point par point » vaut aussi dans n’importe quel groupe multiplicatif cyclique. Pour une dérivation plus algébrique, voir https://www.sciencedirect.com/science/article/pii/S002200007...
    On appelle parfois cela une « FFT harmonique », et il existe aussi des FFT non harmoniques : [LCH14] « additive NTT » sur GF(2^n), [HLP24] circle FFT sur le cercle unité X^2+Y^2=1 d’un corps fini, et [BCKL21] ecfft sur des chaînes d’isogénies de courbes elliptiques
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • Qui a été le premier à proposer d’utiliser la FFT pour une multiplication de polynômes plus rapide ?
    Je me suis posé la question récemment et j’ai cherché ; même si je n’ai pas vraiment réussi à suivre les citations, je suis remonté au moins jusqu’à un article de David Eppstein de 1995 [0]. Il l’utilise pour résoudre efficacement le problème des sommes partielles après des mises à jour incrémentales. Je suis sûr que ça figurait déjà plus tôt dans le TAOCP de Knuth.
    Le fait qu’on puisse aussi résoudre en temps sous-exponentiel le problème exact de la somme de sous-ensemble avec répétitions autorisées grâce à la multiplication de polynômes par FFT m’a pas mal stupéfié [1]. Le point important, c’est que cet algorithme est en O(N log N), où N n’est pas la taille de l’ensemble mais l’élément maximal ; ce n’est donc pas une sorte de contre-exemple à P ≠ NP.
    [0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
    [1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...

  • Je pense que tout apprentissage automatique revient à résoudre des équations de convolution.
    Cet article l’aborde dans le contexte de l’apprentissage par renforcement https://arxiv.org/abs/1712.06115, mais la plupart des approches rentrent dans ce paradigme.

    • En gros, tu veux dire des méthodes à noyau, non ?

  • Je viens d’implémenter un algorithme (matrix profile) qui utilise la FFT pour calculer des produits scalaires sur un grand ensemble de sous-séquences de séries temporelles. La longueur n de la série temporelle peut atteindre plusieurs centaines de millions.
    Le calcul rapide de convolution via FFT fait passer le temps de calcul de O(n) à O(log n), et à cette échelle le gain de vitesse est énorme. Avec un GPU, ça va encore plus vite, par exemple traiter 10 millions de points de données en 0,1 seconde sur un laptop.

  • Le « truc » essentiel de cette opération semble être cette prise de conscience :

    Autrement dit, effectuer la convolution de deux signaux dans le domaine temporel revient à multiplier les deux signaux dans le domaine fréquentiel.
    C’est un bon article, qui découpe une idée complexe en étapes beaucoup plus petites et permet même à quelqu’un de faible en maths comme moi de la comprendre tant bien que mal. Mais est-ce qu’il manque une étape intermédiaire ? Ou bien est-ce laissé au lecteur comme exercice ? À ce stade, j’avais déjà mobilisé tout ce que j’avais comme compétences en maths, et ça m’a un peu fait l’effet de « et maintenant dessinez le reste de ce fichu hibou ». Je suis le seul ? L’article lui-même était vraiment très bon.

    • Je ne sais pas si ça aide, mais : la multiplication de deux polynômes qu’on apprend à l’école est en réalité une convolution.
      Il existe cette propriété selon laquelle « effectuer la convolution de deux signaux dans le domaine temporel revient à multiplier les deux signaux dans le domaine fréquentiel », et la FFT permet de passer du domaine temporel au domaine fréquentiel. On peut donc déplacer les polynômes dans le domaine fréquentiel avec la FFT, puis, dans ce domaine, il suffit de faire une multiplication. C’est plus rapide qu’une convolution. Je me demande si cela rend l’étape manquante plus claire ; s’il manque quelque chose, je peux mettre l’article à jour.

  • Dans ce cas, la factorisation d’entiers est-elle une déconvolution discrète ? Je me demande si, en mettant côte à côte la représentation FFT — c’est-à-dire l’inverse de la multiplication point par point — et tableax, c’est-à-dire la multiplication longue ordinaire / l’addition avec retenues, on obtiendrait assez d’informations, grâce à la rupture de symétrie, pour trouver un algorithme rapide.

  • Bien sûr, la multiplication naïve de polynômes est lente par rapport au degré du polynôme. Mais en pratique, quand a-t-on vraiment besoin de manipuler deux polynômes de degré 100 ?
    C’est pour cela que j’ai l’impression que les systèmes de calcul formel n’utilisent pas ce genre de méthode

    • Les systèmes de calcul formel, par exemple chebfun de Matlab, transforment des fonctions arbitraires en polynômes de degré 100 ou plus afin de trouver plus facilement leurs racines, optimums, etc.
    • C’est très courant dans la correction d’erreurs et le traitement du signal
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... en est un exemple
    • Comme je voulais rétro-ingénierer les paramètres de checksum CRC d’un gros fichier, j’ai créé un programme[1] qui convertit le fichier en un polynôme sur GF(2) de plusieurs millions de degrés et calcule le plus grand commun diviseur. Sans multiplication basée sur la FFT, ce serait impossible dans un délai raisonnable
      [1]: https://github.com/8051enthusiast/delsum
    • Cette perspective par convolution et les kernels GPU rapides pour la FFT ont été utilisés, avant Mamba, dans certains modèles à espace d’états pour la modélisation de longues séquences, où le polynôme est la séquence d’entrée
      Les billets de blog de Hazy Research de 2020 à 2023 contiennent beaucoup d’informations sur cette approche
    • Voir https://news.ycombinator.com/item?id=40306339
      « (…) dans des travaux de recherche en physique, j’ai déjà manipulé des expressions de presque 1 téraoctet, avec plus de 100 millions de termes »