Algorithme de convolution rapide et précis (comme FFT) pour une plage dynamique élevée?

8

Il semble que la convolution basée sur la FFT souffre d'une résolution limitée en virgule flottante en raison de tout évaluer autour des racines de l'unité, comme vous pouvez le voir dans l' erreur de facteur dans ce code Python:1014

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]

Existe-t-il des algorithmes de convolution rapide qui ne souffrent pas de ce problème?
Ou la convolution directe (temps quadratique) est-elle le seul moyen d'obtenir une solution précise?

(Que ces petits nombres soient suffisamment importants pour ne pas être coupés est hors de propos.)

user541686
la source
Notez que convolve()n'appelle que fftconvolve()maintenant, si la taille d'entrée est grande. Spécifiez method='direct'si vous voulez directement.
endolith
@endolith: Bon point! Je viens de l'apprendre récemment, mais je l'ai oublié ici.
user541686

Réponses:

5

Avis de non-responsabilité: je sais que ce sujet est plus ancien, mais si l'on recherche une "plage dynamique élevée de convolution rapide et précise" ou similaire, c'est l'un des premiers des quelques résultats décents. Je veux partager mes idées sur ce sujet afin que cela puisse aider quelqu'un à l'avenir. Je m'excuse si je pourrais utiliser les mauvais termes dans ma réponse, mais tout ce que j'ai trouvé sur ce sujet est assez vague et conduit à la confusion même dans ce fil. J'espère que le lecteur comprendra de toute façon.

La convolution directe est généralement précise à la précision de la machine pour chaque point, c'est-à-dire que l' erreur relative est généralement approximativement ou proche de 1.e-16 pour une double précision pour chaque point du résultat. Chaque point a 16 chiffres corrects. Cependant, les erreurs d'arrondi peuvent être importantes pour des convolutions inhabituellement importantes, et à proprement parler, il faut faire attention à l'annulation et utiliser quelque chose comme la sommation de Kahan et des types de données de précision suffisamment élevée, mais dans la pratique, l'erreur est presque toujours optimale.

L'erreur d'une convolution FFT en dehors des erreurs d'arrondi est une erreur "relative globale", ce qui signifie que l'erreur en chaque point dépend de la précision de la machine et de la valeur de crête du résultat. Par exemple, si la valeur maximale du résultat est 2.e9, alors l'erreur absolue en chaque point est . Donc, si une valeur dans le résultat est supposée être très petite, disons21091016=2107109, l'erreur relative à ce point peut être énorme. La convolution FFT est fondamentalement inutile si vous avez besoin de petites erreurs relatives dans la queue de votre résultat, par exemple vous avez une décroissance quelque peu exponentielle de vos données et avez besoin de valeurs précises dans la queue. Fait intéressant, si la convolution FFT n'est pas limitée par cette erreur, elle a des erreurs d'arrondi beaucoup plus petites que la convolution directe, car vous faites évidemment moins d'ajouts / multiplications. C'est en fait pourquoi les gens prétendent souvent que la convolution FFT est plus précise, et ils ont presque raison dans un certain sens, ils peuvent donc être assez catégoriques.

Malheureusement, il n'y a pas de solution universelle facile pour obtenir des convolutions rapides et précises, mais en fonction de votre problème, il peut y en avoir une ... J'en ai trouvé deux:

Si vous avez des noyaux lisses qui peuvent être bien approchés par un polynôme dans la queue, alors la méthode multipolaire rapide à boîte noire avec interpolation de Chebyshev pourrait être intéressante pour vous. Si votre noyau est "sympa", cela fonctionne parfaitement: vous obtenez à la fois une complexité de calcul linéaire (!) Et une précision de précision de la machine. Si cela correspond à votre problème, vous devez l'utiliser. Ce n'est cependant pas facile à mettre en œuvre.

Pour certains noyaux spécifiques (fonctions convexes je pense, généralement à partir des densités de probabilité), vous pouvez utiliser un "décalage exponentiel" pour obtenir une erreur optimale dans une partie de la queue du résultat. Il y a une thèse de doctorat et un github avec une implémentation de python qui l'utilisent systématiquement, et l'auteur l'appelle une convolution FFT précise . Dans la plupart des cas, cela n'est cependant pas très utile, car soit il régresse en convolution directe, soit vous pouvez utiliser la convolution FFT de toute façon. Bien que le code le fasse automatiquement, ce qui est bien sûr bien sûr.

--------------------ÉDITER:--------------------

J'ai regardé un peu l' algorithme de Karatsuba (j'ai en fait fait une petite implémentation), et il me semble qu'il a généralement un comportement d'erreur similaire à la convolution FFT, c'est-à-dire que vous obtenez une erreur par rapport à la valeur de crête du résultat. En raison de la nature de division et de conquête de l'algorithme, certaines valeurs dans la queue du résultat ont en fait une meilleure erreur, mais je ne vois pas de moyen systématique facile de dire lesquelles ou en tout cas comment utiliser cette observation. Dommage, au début, je pensais que Karatsuba pourrait être quelque chose d'utile entre la convolution directe et la FFT. Mais je ne vois pas de cas d'utilisation courants où Karatsuba devrait être préféré aux deux algorithmes de convolution communs.

Et pour ajouter au décalage exponentiel que j'ai mentionné ci-dessus: Il existe de nombreux cas où vous pouvez l'utiliser pour améliorer le résultat d'une convolution, mais encore une fois, ce n'est pas une solution universelle. J'utilise en fait cela avec la convolution FFT pour obtenir de très bons résultats (dans le cas général pour toutes les entrées: au pire, même erreur que la convolution FFT normale, au mieux erreur relative à chaque point de la précision de la machine). Mais encore une fois, cela ne fonctionne vraiment bien que pour des noyaux et des données spécifiques, mais pour moi, le noyau et les données ou quelque peu exponentiels en décomposition.

oli
la source
+1 bienvenue et merci beaucoup d'avoir posté ceci! :)
user541686
1
Hou la la! j'ai aussi appris quelque chose et c'est un nouveau terme pour quelque chose que je fais depuis 1993. cet algorithme de sommation de Kahan semble être exactement le même que ce que j'appelais la mise en forme du bruit avec un zéro dans la fonction de transfert de bruit à la sortie placé juste à DC ou le zéro est placé à sur le plan . Randy Yates l'a appelé " économie de fraction ", qui est un nom générique concis. je me demande qui est mr / ms Kahan et quand cela est crédité. z=1z
robert bristow-johnson
2
La publication originale de Kahan semble être de 1964.
oli
c'est la surprise d'aujourd'hui. En fait, un certain temps @DanBoschen avait demandé un puzzle dsp, compte tenu de la plage dynamique des nombres à virgule flottante, qui concernait en fait ce même concept d'ajouter des très petits nombres aux très grands nombres ...
Fat32
3

Un des candidats est l' algorithme de Karatsuba , qui s'exécute en temps . Ce n'est pas basé sur la transformation. Il y a aussi du code avec une explication dans les archives de code source de Music-DSP, qui ressemble à une découverte indépendante d'un algorithme similaire.O(Nlog23)O(N1.5849625)

Le test d'une implémentation Python de l'algorithme Karatsuba (installé par sudo pip install karatsuba) à l'aide des nombres dans votre question montre que même avec des nombres à virgule flottante 64 bits, l'erreur relative est importante pour l'une des valeurs de sortie:

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])

qui imprime:

Karatsuba:
[1.0, 1.9984014443252818e-15, 1.0000000000000001e-30]
Direct:
[1.00000000000000000E+00 2.00000000000000016E-15 1.00000000000000008E-30]
Olli Niemitalo
la source
2
Il y a un supplément] à la fin du lien vers l' algorithme de Karatsuba
+1 parce que c'est génial et qu'il ne m'était jamais venu à l'esprit que Karatsuba était un algorithme de convolution, mais ce serait bien si vous pouviez expliquer pourquoi cela devrait résoudre ce problème. Je peux facilement le voir pour le cas 2x2, mais dans le paramètre récursif général, je ne vois pas pourquoi cela devrait résoudre ce problème. Il me semble plausible que ce ne soit même pas réparable en général, mais je ne sais pas.
user541686
1
@OlliNiemitalo: Eh bien, la façon la plus simple de l'expliquer est que je veux que l'erreur relative soit faible par rapport à la convolution directe . (Toute définition raisonnable de "faible" fonctionnerait ici ... l'erreur relative que je reçois avec FFT est comme qui n'est pas faible par définition.)O(n2)1014
user541686
1
Les doubles IEEE n'ont qu'une précision de 15 à 16 chiffres décimaux dans le cas général. Donc 1e-14 est une erreur de taille raisonnable pour une séquence d'un certain nombre d'opérations arithmétiques (sauf si vous choisissez quelques valeurs magiques).
hotpaw2
1
Si vous avez déjà conçu un additionneur à virgule flottante, vous saurez que l'exposant est déterminé par le résultat de la mantisse pendant la normalisation. Vous avez choisi des nombres qui produisent une mantisse étroite improbable.
hotpaw2
3

Plutôt que de supprimer l'algorithme de convolution rapide, pourquoi ne pas utiliser une FFT avec une plage dynamique plus élevée?

Une réponse à cette question montre comment utiliser la bibliothèque Eigen FFT avec boost multiprecision.

Mark Borgerding
la source
2

Je crois que la précision de l'algorithme Cordic peut être étendue autant que vous le souhaitez, si c'est le cas, utilisez un DFT entier et une longueur de mot appropriée à votre problème.

La même chose serait vraie avec une convolution directe, utilisez des entiers très longs.


la source
1

La convolution temporelle quadratique pour obtenir un résultat DFT est généralement moins précise (peut générer un bruit numérique de quantification plus fin, en raison d'une superposition plus profonde des étapes arithmétiques) que l'algorithme FFT typique lors de l'utilisation des mêmes types arithmétiques et unités de fonctionnement.

Vous voudrez peut-être essayer des types de données de précision plus élevée (précision quadruple ou arithmétique bignum).

hotpaw2
la source
Er, c'est en utilisant les mêmes types arithmétiques et les unités d'exploitation est pas? C'est clairement plus précis. Je pense que le type de bruit dont vous parlez n'est pas le même que celui dont je parle. Les racines de l'unité ont une magnitude de 1, ce qui signifie qu'elles ne peuvent tout simplement pas représenter de très petites valeurs. Cela ne semble pas totalement lié à la question de savoir comment le bruit se propage à travers le système.
user541686
Cela ne semble plus précis dans votre exemple, car vous avez choisi une longueur et des valeurs où l'arrondi s'est avéré être en votre faveur. Essayez une gamme de convolutions beaucoup plus longues avec beaucoup plus de coefficients non nuls avec une distribution contenant un large ordre de grandeur.
hotpaw2
Le problème que j'essaie de résoudre n'a cependant rien à voir avec l'arrondi. C'est un problème différent que je n'essaie pas de résoudre. Les exemples originaux que j'avais étaient exactement comme ce que vous venez de dire, et ils fonctionnaient très bien avec convolution directe mais ont été détruits par FFT.
user541686
L'arrondi (ou d'autres méthodes de quantification) est impliqué dans toute l'arithmétique de précision finie. Certains résultats de calcul changent lorsqu'ils sont arrondis, d'autres non ou changent moins.
hotpaw2
Je n'ai jamais prétendu le contraire. Ce que je viens de vous dire, c'est que le problème que j'essaie de résoudre n'a rien à voir avec l'arrondissement. C'est un problème différent. Je ne me soucie pas d'éviter l'arrondi, mais je me soucie d'éviter ce problème.
user541686