Faltung, schnelle Fourier-Transformation und Polynome (2022)

 (alvarorevuelta.com)
1 Punkte von GN⁺ 2024-07-02 | 1 Kommentare | Auf WhatsApp teilen
  • Wenn man Polynome hohen Grades nach der Schulmethode ausmultipliziert, müssen alle Paare von Termen multipliziert werden; die Kosten von O(n²) werden daher schnell zum Flaschenhals
  • Die Multiplikation von Koeffizientenvektoren eines Polynoms entspricht der Faltung diskreter Signale; das Ergebnis von [2, 3, 4] und [5, 6, 7] ist [10, 27, 52, 45, 28]
  • Die DFT überführt ein diskretes Signal in den Frequenzbereich, und die FFT berechnet dieselbe Transformation in O(n log n), was bei großen Eingaben den Unterschied macht
  • Eine Faltung im Zeitbereich wird im Frequenzbereich zu einer elementweisen Multiplikation. Wandelt man also per FFT um, multipliziert und wandelt anschließend per IFFT zurück, lässt sich die Polynommultiplikation schneller durchführen
  • Bei kleinen Graden können die Kosten für Hin- und Rücktransformation per FFT/IFFT den Vorteil aufheben, doch mit wachsendem Grad wird das FFT-Verfahren effizienter

Warum Polynommultiplikation langsam wird

  • Ein Polynom P(x) wird als Summe von Koeffizienten a_k und Potenzen der Variablen x dargestellt
    • Das Beispiel P(x)=5x²+2x+9 ist ein Polynom zweiten Grades
    • Der Koeffizientenvektor kann je nach Schreibweise etwa als [5, 2, 9] oder [9, 2, 5] dargestellt werden
  • Addition und Subtraktion sind vergleichsweise einfach, weil man nur Terme gleichen Grades addiert oder subtrahiert
    • In Python kann man mit zip(p, q) die einzelnen Koeffizienten durchlaufen und a + b oder a - b berechnen
    • Bei unterschiedlichen Graden kann zip_longest verwendet werden
  • Bei der Multiplikation muss jeder Term mit jedem anderen multipliziert und anschließend müssen Terme gleichen Grades wieder zusammengefasst werden, wodurch der Rechenaufwand steigt
    • Das Ergebnis von (2x²+3x+4) × (5x²+6x+7) ist 10x⁴+27x³+52x²+45x+28
    • Die Komplexität dieses Verfahrens beträgt O(n²), und mit wachsendem Grad steigt die Zahl der benötigten Multiplikationen

Koeffizientenvektoren und Faltung

  • Im diskreten Bereich ist die Faltung zweier Signale p und q definiert als y[n]=Σ p[k]·q[n-k]
  • Die Berechnung erfolgt, indem q umgedreht und von links nach rechts über p verschoben wird, wobei die Produkte der überlappenden Elemente addiert werden
  • Die Beispielsignale sind wie folgt
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • Wenn man q umdreht und verschiebt, entstehen die einzelnen Ausgabekoeffizienten in dieser Reihenfolge
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • Das Faltungsergebnis ist y = [10, 27, 52, 45, 28]
    • Es entspricht den Koeffizienten von 10x⁴+27x³+52x²+45x+28, die man durch Polynommultiplikation erhält
    • Daher lässt sich Polynommultiplikation als Faltung von Koeffizientenvektoren auffassen

Fourier-Transformation und FFT

  • Die Fourier-Transformation wandelt ein Signal vom Zeitbereich in den Frequenzbereich um
    • Aus zeitlicher Sicht betrachtet man ein Signal über seine Werte zu bestimmten Zeitpunkten
    • Aus Frequenzsicht interpretiert man ein Signal als Summe verschiedener Schwingungsfrequenzen
  • Schwingungsfrequenzen werden durch Sinus und Kosinus dargestellt und haben jeweils Koeffizient und Phase
  • Wendet man die FFT auf eine reine 5-Hz-Sinuswelle an, erscheint sie im Frequenzbereich deltaartig an der 5-Hz-Position
    • Das zeigt, dass die Sinuswelle im Zeitbereich durch einen einzelnen 5-Hz-Sinus dargestellt werden kann
  • Die einschlägigen Begriffe werden wie folgt unterschieden
    • Fourier Transform(FT): im kontinuierlichen Bereich definierte Fourier-Transformation
    • Discrete Fourier Transform(DFT): für diskrete Signale definierte Fourier-Transformation
    • Fast Fourier Transform(FFT): Algorithmus, der die DFT statt in O(n²) in O(n log n) berechnet
  • Die DFT wandelt ein diskretes Zeitsignal x[n] in X[k] im Frequenzbereich um
    • Jedes X[k] wird berechnet, indem die Eingabesamples mit einer komplexen Zahl multipliziert und aufsummiert werden, die eine bestimmte Frequenz repräsentiert

Umwandlung in Multiplikation im Frequenzbereich

  • Der zentrale Vorteil der DFT und des Frequenzbereichs besteht darin, dass sich Faltung in elementweise Multiplikation umwandeln lässt
    • Zwei Signale im Zeitbereich zu falten ist gleichbedeutend damit, die beiden Signale im Frequenzbereich zu multiplizieren
    • Multiplikation lässt sich schneller berechnen als Faltung
  • Das Verfahren zur schnellen Polynommultiplikation sieht wie folgt aus
    • Die Polynome per FFT in den Frequenzbereich transformieren: O(n log n)
    • Im Frequenzbereich elementweise multiplizieren: O(n)
    • Das Ergebnis per IFFT wieder in den Zeitbereich zurücktransformieren: O(n log n)
  • Insgesamt kann Polynommultiplikation mit der FFT in O(n log n) durchgeführt werden
  • Bei großen Polynomen ist das schneller als die Schulmethode mit O(n²)

Python-Implementierung und Benchmark

  • multiply_naive multipliziert mit einer doppelten Schleife alle Koeffizientenpaare und addiert sie an der Ergebnisposition i + j
    • Die Ergebnislänge ist len(p) + len(q) - 1
    • Die Komplexität beträgt O(n²)
  • multiply_fft führt die Koeffizientenmultiplikation auf Basis von FFT/IFFT aus
    • Es berechnet eine Zweierpotenz-Länge, die mindestens len(p) + len(q) - 1 aufnehmen kann
    • Die beiden Eingaben werden mit np.pad gepaddet
    • Die mit np.fft.fft transformierten Werte werden elementweise multipliziert
    • Nach der Rücktransformation mit np.fft.ifft wird der Realteil gerundet und in ganzzahlige Koeffizienten umgewandelt
  • Für die Beispieleingaben p = [2, 3, 4], q = [5, 6, 7] liefern beide Verfahren [10, 27, 52, 45, 28]
  • Im Benchmark wird multiply_convolve, das np.convolve verwendet, mit dem FFT-Verfahren verglichen, statt multiply_naive zu nutzen
    • Denn multiply_naive ist wegen der Python-Schleifen langsam und daher schwer direkt mit dem FFT-Verfahren zu vergleichen, das np.fft.fft verwendet
    • np.convolve führt dieselbe Operation in Low-Level-C-Code aus
  • Der Grad wird im Bereich range(1, 30000, 1000) erhöht, und für jeden Grad werden zwei Polynome mit zufälligen Koeffizienten zwischen 1 und 999999 erzeugt
    • Für jedes Verfahren wird mit n_runs = 5 die Durchschnittszeit gemessen
    • Bei niedrigen Graden kann das FFT-Verfahren wegen der Kosten für Hin- und Rücktransformation per FFT/IFFT im Nachteil sein
    • Mit steigendem Grad zeigt das FFT-Verfahren deutlich effizientere Ergebnisse

Verwandte Beiträge

1 Kommentare

 
GN⁺ 2024-07-02
Kommentare auf Hacker News

  • Was mich an solchen Erklärungen immer stört, ist, dass sie meist numerische Fehler vergessen.
    Koeffizientenmultiplikation kann man nicht einfach als „konstante Zeit“ abstrahieren. Wenn man das tut, könnte man genauso gut gleich die gesamte Multiplikation abstrahieren. Berücksichtigt man die numerische Genauigkeit, liegt es eher bei O(n (log n)^3) [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • Die in diesem Artikel zitierte Fehlerschranke ist übermäßig pessimistisch. In der neuesten Ausgabe von Knuth steht die korrekte Schranke, weil ich ihn darauf hingewiesen habe.
    • Es wäre schön, wenn man die im OP-Artikel erwähnten quaternionenbasierten Operationen nutzen könnte, um Multiplikationsfehler zu reduzieren oder sogar ganz zu eliminieren [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
    • Wenn die Koeffizienten ganzzahlig sind, kann man mit einer NTT mit ausreichend großem Modulus exakte Ergebnisse erhalten; insbesondere in Hardware kann auch die Multiplikationszeit kürzer ausfallen.
    • Deshalb unterscheidet man zwischen Informatik und Software Engineering :)

  • Mit dieser Methode kann man lange Zahlen miteinander multiplizieren. Der Kern ist, dass Polynommultiplikation dasselbe ist wie die normale schriftliche Multiplikation langer Zahlen ohne Überträge (carry).
    Wenn man zum Beispiel eine 1000-stellige Zahl hat, nimmt man jede Ziffer als Koeffizienten eines Polynoms mit 1000 Elementen. Anschließend kann man diese Polynome mit der im Artikel beschriebenen FFT-Methode multiplizieren. Um das Ergebnis wieder in eine Zahl umzuwandeln, muss man die Überträge verarbeiten. Wenn ein Element größer als 10 ist, gibt man den Überschuss an die nächste Stelle weiter und wandelt die Koeffizienten in Ziffern um.
    Das ist die Grundidee; bei der für die Überträge nötigen Genauigkeit und der Garantie, dass das Runden der FFT-Ergebnisse auf die nächste ganze Zahl korrekt ist, gibt es einige Feinheiten. Auf diese Weise führt GMP, die führende Bibliothek in diesem Bereich, Multiplikationen großer Zahlen durch.

    • Wie gesagt, Dezimalzahlen lassen sich als Polynom mit x=10 darstellen, das ist nachvollziehbar. Zum Beispiel ist 983 = 9x^2 + 8x + 3, also [9, 8, 3].
      Ich frage mich, wie groß Zahlen sein müssen, damit das in der Praxis sinnvoll ist, und wofür es eingesetzt wird.

  • Falls du es noch nicht gesehen hast, ist dieses Video empfehlenswert:
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    Es leitet den FFT-Algorithmus aus der Polynommultiplikation her und ist wirklich hervorragend. Ich schaue es mir etwa alle sechs Monate wieder an.

  • Die FFT-Eigenschaft „Faltung ist punktweise Multiplikation“ gilt auch für beliebige zyklische multiplikative Gruppen. Für eine stärker algebraische Herleitung siehe https://www.sciencedirect.com/science/article/pii/S002200007...
    Das wird manchmal als „harmonische FFT“ bezeichnet; es gibt auch nicht-harmonische FFTs: [LCH14] „additive NTT“ über GF(2^n), [HLP24] circle FFT auf dem Einheitskreis X^2+Y^2=1 eines endlichen Körpers, [BCKL21] ecfft über Isogeniefolgen elliptischer Kurven.
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • Wer hat eigentlich zuerst vorgeschlagen, FFT für schnellere Polynom-Multiplikation zu verwenden?
    Ich habe kürzlich aus Neugier recherchiert und bin, auch wenn ich die Zitationskette nicht sauber verfolgen konnte, bis zu einem Paper von David Eppstein aus dem Jahr 1995 [0] zurückgegangen. Dort wird sie verwendet, um nach inkrementellen Aktualisierungen das Partialsummenproblem effizient zu lösen. In Knuths TAOCP stand das bestimmt schon früher
    Dass man mit FFT-Polynom-Multiplikation das exakte Partialsummenproblem mit Wiederholungen in subexponentieller Zeit lösen kann, fand ich ebenfalls ziemlich verblüffend [1]. Wichtig ist, dass dieser Algorithmus O(N log N) ist, wobei N nicht die Größe der Menge, sondern das maximale Element ist; es ist also kein Gegenbeispiel zu P ≠ NP
    [0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
    [1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...

  • Ich denke, dass alles maschinelle Lernen letztlich darin besteht, Faltungsgleichungen zu lösen
    Dieses Paper behandelt es im Kontext von Reinforcement Learning https://arxiv.org/abs/1712.06115, aber die meisten Ansätze passen in dieses Paradigma

    • Heißt das im Grunde nicht Kernel-Methoden?

  • Ich habe gerade einen Algorithmus (matrix profile) implementiert, der FFT verwendet, um Skalarprodukte für eine große Menge von Teilsequenzen einer Zeitreihe zu berechnen. Die Länge n der Zeitreihe kann bis in den Bereich von Hunderten Millionen gehen
    Durch schnelle Faltungsberechnung per FFT sinkt die Rechenzeit von O(n) auf O(log n), und in dieser Größenordnung ist der Geschwindigkeitsgewinn enorm. Mit GPU wird es noch schneller, etwa 10 Millionen Datenpunkte in 0,1 Sekunden auf einem Laptop

  • Der zentrale „Trick“ dieser Operation scheint diese Erkenntnis zu sein:

    Mit anderen Worten: Eine Faltung zweier Signale im Zeitbereich durchzuführen ist dasselbe, wie die beiden Signale im Frequenzbereich zu multiplizieren.
    Ein guter Artikel, der eine komplexe Idee in viel kleinere Schritte zerlegt und sie dadurch sogar für jemanden wie mich, der in Mathe schwach ist, irgendwie verständlich macht. Aber wurde da ein Zwischenschritt ausgelassen? Oder als Übungsaufgabe für die Lesenden stehen gelassen? Bis zu diesem Punkt hatte ich meine mathematischen Fähigkeiten schon maximal ausgereizt, und es fühlte sich ein bisschen an wie „und jetzt zeichne den restlichen verdammten Vogel“. Geht es nur mir so? Der Artikel selbst war wirklich gut

    • Ob es hilft, weiß ich nicht, aber: Die Multiplikation zweier Polynome, wie man sie in der Schule lernt, ist tatsächlich eine Faltung
      Es gibt die Eigenschaft: „Eine Faltung zweier Signale im Zeitbereich durchzuführen ist dasselbe, wie die beiden Signale im Frequenzbereich zu multiplizieren“, und FFT ermöglicht es, vom Zeitbereich in den Frequenzbereich zu transformieren. Deshalb verschiebt man die Polynome per FFT in den Frequenzbereich und muss dort nur noch multiplizieren. Das ist schneller als eine Faltung. Ich frage mich, ob der fehlende Schritt damit klarer wird; falls noch etwas fehlt, kann der Artikel aktualisiert werden

  • Ist Ganzzahlfaktorisierung dann diskrete Dekonvolution? Ich frage mich, ob man genug Informationen für einen schnellen Algorithmus gewinnen könnte, wenn man die FFT-Darstellung, also die Umkehrung der punktweisen Multiplikation, und tableax, also die normale schriftliche Multiplikation/Addition mit Übertrag, nebeneinanderstellt und dadurch die Symmetrie bricht

  • Natürlich ist naive Polynommultiplikation in Bezug auf den Grad des Polynoms langsam. Aber wann muss man in der Praxis schon einmal mit zwei Polynomen 100. Grades umgehen?
    Aus diesem Grund hat man den Eindruck, dass Computer-Algebra-Systeme diese Methode nicht verwenden

    • Computer-Algebra-Systeme, etwa Matlab chebfun, wandeln beliebige Funktionen in Polynome 100. Grades oder höher um, um Nullstellen, Optima usw. leichter zu finden
    • Bei Fehlerkorrektur und Signalverarbeitung ist das sehr verbreitet
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... ist ein Beispiel dafür
    • Weil ich die CRC-Prüfsummenparameter großer Dateien rückentwickeln wollte, habe ich ein Programm[1] geschrieben, das Dateien in GF(2)-Polynome mit Millionen von Graden umwandelt und den größten gemeinsamen Teiler berechnet. Ohne FFT-basierte Multiplikation ist das in angemessener Zeit nicht möglich
      [1]: https://github.com/8051enthusiast/delsum
    • Diese Faltungs-Perspektive und schnelle GPU-Kernels für FFT wurden vor Mamba in einigen Zustandsraummodellen zur Modellierung langer Sequenzen verwendet; dabei sind die Polynome die Eingabesequenzen
      In den Blogbeiträgen von Hazy Research aus den Jahren 2020 bis 2023 gibt es viele Informationen zu diesem Ansatz
    • Siehe https://news.ycombinator.com/item?id=40306339
      „(...) In der Physikforschung habe ich einmal mit Ausdrücken gearbeitet, die fast 1 Terabyte groß waren und mehr als 100 Millionen Terme hatten“