Schneller und genauer Faltungsalgorithmus (wie FFT) für hohen Dynamikbereich?

Es scheint, dass die FFT-basierte Faltung unter einer begrenzten Gleitkommaauflösung leidet, da alles um die Wurzeln der Einheit herum ausgewertet wird, wie Sie im -Faktorfehler in diesem Python-Code sehen können:$10^{14}$

from scipy.signal import convolve, fftconvolve
a = [1.0, 1E-15]
b = [1.0, 1E-15]
convolve(a, b)     # [  1.00000000e+00,   2.00000000e-15,   1.00000000e-30]
fftconvolve(a, b)  # [  1.00000000e+00,   2.11022302e-15,   1.10223025e-16]

Gibt es schnelle Faltungsalgorithmen, die nicht unter diesem Problem leiden?
Oder ist eine direkte (quadratische) Faltung der einzige Weg, um eine genaue Lösung zu erhalten?

(Ob solch kleine Zahlen signifikant genug sind, um nicht abzuhacken, ist neben meinem Standpunkt.)

Haftungsausschluss: Ich weiß, dass dieses Thema älter ist, aber wenn man nach einem "schnellen, genauen Faltungs-Hochdynamikbereich" oder ähnlichem sucht, ist dies eines der ersten von nur wenigen anständigen Ergebnissen. Ich möchte meine Erkenntnisse zu diesem Thema teilen, damit es in Zukunft jemandem helfen kann. Ich entschuldige mich, wenn ich in meiner Antwort möglicherweise die falschen Begriffe verwende, aber alles, was ich zu diesem Thema gefunden habe, ist ziemlich vage und führt selbst in diesem Thread zu Verwirrung. Ich hoffe der Leser wird es trotzdem verstehen.

Die direkte Faltung ist meistens auf die Maschinengenauigkeit für jeden Punkt genau, dh der relative Fehler liegt normalerweise ungefähr oder nahe bei 1.e-16 für die doppelte Genauigkeit für jeden Punkt des Ergebnisses. Jeder Punkt hat 16 korrekte Ziffern. Rundungsfehler können für untypisch große Windungen von Bedeutung sein, und genau genommen sollte man beim Löschen vorsichtig sein und so etwas wie Kahan-Summierung und ausreichend genaue Datentypen verwenden, aber in der Praxis ist der Fehler fast immer optimal.

2.e9$2\cdot10^9\cdot10^{-16} = 2\cdot10^{-7}$$10^{-9}$kann der relative Fehler in diesem Punkt sehr groß sein. Die FFT-Faltung ist grundsätzlich nutzlos, wenn Sie kleine relative Fehler im Ende Ihres Ergebnisses benötigen, z. B. wenn Ihre Daten etwas exponentiell abfallen und genaue Werte im Ende benötigen. Interessanterweise weist die FFT-Faltung, wenn sie nicht durch diesen Fehler begrenzt ist, im Vergleich zur direkten Faltung viel kleinere Rundungsfehler auf, da Sie offensichtlich weniger Additionen / Multiplikationen durchführen. Dies ist tatsächlich der Grund, warum Menschen oft behaupten, dass die FFT-Faltung genauer ist, und sie haben in gewissem Sinne fast recht, so dass sie ziemlich unnachgiebig sein können.

Leider gibt es keine einfache universelle Lösung , um schnelle und genaue Windungen zu erhalten, aber abhängig von Ihrem Problem kann es eine geben ... Ich habe zwei gefunden:

Wenn Sie glatte Kernel haben, die durch ein Polynom im Schwanz gut angenähert werden können, ist die Black-Box-Fast-Multipole-Methode mit Chebyshev-Interpolation möglicherweise für Sie interessant. Wenn Ihr Kernel "nett" ist, funktioniert dies tatsächlich perfekt: Sie erhalten sowohl lineare (!) Rechenkomplexität als auch Maschinengenauigkeit. Wenn dies zu Ihrem Problem passt, sollten Sie es verwenden. Es ist jedoch nicht einfach zu implementieren.

Für einige spezifische Kernel (konvexe Funktionen, glaube ich, normalerweise aus Wahrscheinlichkeitsdichten) können Sie eine "exponentielle Verschiebung" verwenden, um einen optimalen Fehler in einem Teil des Endes des Ergebnisses zu erhalten. Es gibt eine Doktorarbeit und einen Github mit einer Python-Implementierung , die diese systematisch verwendet, und der Autor nennt sie eine genaue FFT-Faltung . In den meisten Fällen ist dies jedoch nicht besonders nützlich, da es entweder zur direkten Faltung zurückkehrt oder Sie die FFT-Faltung trotzdem verwenden können. Obwohl der Code es automatisch macht, ist das natürlich schön.

--------------------BEARBEITEN:--------------------

Ich habe mir den Karatsuba- Algorithmus ein wenig angesehen (ich habe tatsächlich eine kleine Implementierung vorgenommen), und für mich sieht es so aus, als ob er normalerweise ein ähnliches Fehlerverhalten wie die FFT-Faltung aufweist, dh Sie erhalten einen Fehler relativ zum Spitzenwert des Ergebnisses. Aufgrund der Teilung und Eroberung des Algorithmus weisen einige Werte am Ende des Ergebnisses tatsächlich bessere Fehler auf, aber ich sehe keine einfache systematische Methode, um festzustellen, welche oder auf jeden Fall, wie diese Beobachtung verwendet werden soll. Schade, zuerst dachte ich, Karatsuba könnte etwas Nützliches zwischen direkter und FFT-Faltung sein. Ich sehe jedoch keine häufigen Anwendungsfälle, in denen Karatsuba den beiden üblichen Faltungsalgorithmen vorgezogen werden sollte.

Und um die oben erwähnte exponentielle Verschiebung zu ergänzen : Es gibt viele Fälle, in denen Sie damit das Ergebnis einer Faltung verbessern können, aber es ist wiederum keine universelle Lösung. Ich verwende dies tatsächlich zusammen mit der FFT-Faltung, um ziemlich gute Ergebnisse zu erzielen (im allgemeinen Fall für alle Eingaben: im schlimmsten Fall der gleiche Fehler wie bei der normalen FFT-Faltung, bestenfalls relativer Fehler in jedem Punkt zur Maschinengenauigkeit). Aber auch dies funktioniert nur für bestimmte Kernel und Daten wirklich gut, aber für mich sowohl Kernel als auch Daten oder etwas exponentiell im Zerfall.

$O\big(N^{\log_23}\big) \approx O\big(N^{1.5849625}\big)$

Das Testen einer Python-Implementierung des Karatsuba-Algorithmus (installiert von sudo pip install karatsuba) unter Verwendung der Zahlen in Ihrer Frage zeigt, dass selbst bei 64-Bit-Gleitkommazahlen der relative Fehler für einen der Ausgabewerte groß ist:

import numpy as np
from karatsuba import *
k = make_plan(range(2), range(2))
l = [np.float64(1), np.float64(1E-15)]
np.set_printoptions(formatter={'float': lambda x: format(x, '.17E')})
print "Karatsuba:"
print(k(l, l)[0:3])
print "Direct:"
print(np.convolve(l, l)[0:3])

welche druckt:

Karatsuba:
[1.0, 1.9984014443252818e-15, 1.0000000000000001e-30]
Direct:
[1.00000000000000000E+00 2.00000000000000016E-15 1.00000000000000008E-30]

Warum nicht eine FFT mit einem höheren Dynamikbereich verwenden, anstatt den Algorithmus für schnelle Faltung zu verschrotten?

Eine Antwort auf diese Frage zeigt, wie die Eigen-FFT-Bibliothek mit Boost-Multipräzision verwendet wird.

Ich glaube, dass die Genauigkeit des Cordic-Algorithmus so weit erweitert werden kann, wie Sie möchten, wenn Sie eine ganzzahlige DFT und eine Wortlänge verwenden, die Ihrem Problem entspricht.

Das gleiche gilt für die direkte Faltung, verwenden Sie sehr lange ganze Zahlen.

Die quadratische Zeitfaltung zum Erhalten eines DFT-Ergebnisses ist normalerweise weniger genau (kann aufgrund einer tieferen Schichtung von arithmetischen Schritten ein endlicheres numerisches Quantisierungsrauschen verursachen) als der typische FFT-Algorithmus, wenn dieselben arithmetischen Typen und Operationseinheiten verwendet werden.

Möglicherweise möchten Sie Datentypen mit höherer Genauigkeit (Quad-Genauigkeit oder Bignum-Arithmetik) ausprobieren.