Algoritmo de convolución rápido y preciso (como FFT) para un alto rango dinámico?

8

Parece que la convolución basada en FFT adolece de una resolución de punto flotante limitada debido a la evaluación de todo en torno a las raíces de la unidad, como se puede ver en el 1014-factor de error en este código de Python:

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]

¿Hay algún algoritmo de convolución rápida que no sufra este problema?
¿O es la convolución directa (tiempo cuadrático) la única forma de obtener una solución precisa?

(El hecho de que tales números pequeños sean lo suficientemente significativos como para no cortarlos está fuera de mi punto).

usuario541686
fuente
Tenga en cuenta que convolve()solo llama fftconvolve()ahora, si el tamaño de entrada es grande. Especifica method='direct'si quieres directo.
Endolith
@endolith: ¡Buen punto! Acabo de enterarme de eso recientemente, pero lo olvidé aquí.
user541686

Respuestas:

5

Descargo de responsabilidad: Sé que este tema es más antiguo, pero si uno está buscando un "rango dinámico alto de convolución rápido y preciso" o similar, este es uno de los primeros de unos pocos resultados decentes. Quiero compartir mis ideas sobre este tema para que pueda ayudar a alguien en el futuro. Pido disculpas si podría usar los términos incorrectos en mi respuesta, pero todo lo que encontré sobre este tema es bastante vago y puede llevar a confusión incluso en este hilo. Espero que el lector lo entienda de todos modos.

La convolución directa es en su mayoría precisa para la precisión de la máquina para cada punto, es decir, el error relativo suele ser aproximadamente o cercano a 1.e-16 para obtener una precisión doble para cada punto del resultado. Cada punto tiene 16 dígitos correctos. Sin embargo, los errores de redondeo pueden ser significativos para convoluciones anormalmente grandes, y estrictamente hablando, uno debe tener cuidado con la cancelación y usar algo como la suma de Kahan y los tipos de datos de precisión suficientemente alta, pero en la práctica el error es casi siempre óptimo.

El error de una convolución FFT aparte de los errores de redondeo es un error "relativo global", lo que significa que el error en cada punto depende de la precisión de la máquina y del valor máximo del resultado. Por ejemplo, si el valor máximo del resultado es 2.e9, entonces el error absoluto en cada punto es21091016=2107. Entonces, si se supone que un valor en el resultado es muy pequeño, digamos109, el error relativo en ese punto puede ser enorme. La convolución FFT es básicamente inútil si necesita pequeños errores relativos en la cola de su resultado, por ejemplo, tiene una disminución algo exponencial de sus datos y necesita valores precisos en la cola. Curiosamente, si la convolución FFT no está limitada por ese error, tiene errores de redondeo mucho más pequeños en comparación con la convolución directa, ya que obviamente hace menos adiciones / multiplicaciones. Esta es realmente la razón por la cual las personas a menudo afirman que la convolución FFT es más precisa, y casi tienen razón en algún sentido, por lo que pueden ser bastante inflexibles.

Desafortunadamente, no hay una solución universal fácil para obtener convoluciones rápidas y precisas, pero dependiendo de su problema puede haber una ... He encontrado dos:

Si tiene granos lisos que se pueden aproximar bien por un polinomio en la cola, entonces el Método rápido multipolar de caja negra con interpolación de Chebyshev podría ser interesante para usted. Si su núcleo es "agradable", esto funciona realmente perfectamente: obtiene tanto la complejidad computacional lineal (!) Como la precisión de la precisión de la máquina. Si esto encaja con su problema, debe usarlo. Sin embargo, no es fácil de implementar.

Para algunos núcleos específicos (funciones convexas, creo, generalmente de densidades de probabilidad) puede usar un "cambio exponencial" para obtener un error óptimo en alguna parte de la cola del resultado. Hay una tesis de doctorado y un github con una implementación de Python que lo usa sistemáticamente, y el autor lo llama convolución FFT precisa . Sin embargo, en la mayoría de los casos esto no es súper útil, ya que regresa a la convolución directa o puede usar la convolución FFT de todos modos. Aunque el código lo hace automáticamente, lo cual es bueno, por supuesto.

--------------------EDITAR:--------------------

Miré un poco el algoritmo de Karatsuba (en realidad hice una pequeña implementación), y para mí parece que generalmente tiene un comportamiento de error similar al de la convolución FFT, es decir, obtienes un error relativo al valor máximo del resultado. Debido a la naturaleza de dividir y conquistar del algoritmo, algunos valores en la cola del resultado en realidad tienen un mejor error, pero no veo una manera sistemática fácil de decir cuáles o, en cualquier caso, cómo usar esta observación. Lástima, al principio pensé que Karatsuba podría ser algo útil entre convolución directa y FFT. Pero no veo casos de uso comunes donde Karatsuba deba preferirse sobre los dos algoritmos de convolución comunes.

Y para agregar al cambio exponencial que mencioné anteriormente: hay muchos casos en los que puede usarlo para mejorar el resultado de una convolución, pero nuevamente no es una solución universal. De hecho, utilizo esto junto con la convolución FFT para obtener resultados bastante buenos (en el caso general de todas las entradas: en el peor error que la convolución FFT normal, en el mejor error relativo en cada punto para la precisión de la máquina). Pero, de nuevo, esto solo funciona muy bien para núcleos y datos específicos, pero para mí tanto el núcleo como los datos o algo exponencial en decadencia.

oli
fuente
¡+1 bienvenido y muchas gracias por publicar esto! :)
user541686
1
¡Guauu! también aprendí algo y ese es un nuevo término para algo que he estado haciendo desde 1993. Este algoritmo de suma de Kahan parece ser exactamente el mismo que había estado llamando forma de ruido con un cero en la función de transferencia de ruido a salida colocado justo en DC o el cero se coloca enz=1 sobre el zavión. Randy Yates lo llamó " ahorro de fracción ", que es un nombre genérico conciso. Me pregunto quién es el señor / a Kahan y cuándo se acredita esto.
Robert Bristow-Johnson
2
La publicación original de Kahan parece ser de 1964.
oli
Es la sorpresa de hoy. En realidad, hace un tiempo, @DanBoschen había pedido un rompecabezas dsp, considerando el rango dinámico de los números de coma flotante, que en realidad se trataba del mismo concepto de agregar números muy pequeños a números muy grandes ...
Fat32
3

Un candidato es el algoritmo Karatsuba , que se ejecuta enO(Nlog23)O(N1.5849625)hora. No está basado en la transformación. También hay un código con una explicación en el Archivo de código fuente Music-DSP, que parece un descubrimiento independiente de un algoritmo similar.

Probar una implementación de Python del algoritmo Karatsuba (instalado por sudo pip install karatsuba) usando los números en su pregunta muestra que incluso con números de coma flotante de 64 bits, el error relativo es grande para uno de los valores de salida:

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

que imprime:

Karatsuba:
[1.0, 1.9984014443252818e-15, 1.0000000000000001e-30]
Direct:
[1.00000000000000000E+00 2.00000000000000016E-15 1.00000000000000008E-30]
Olli Niemitalo
fuente
2
Hay un extra] al final del enlace al algoritmo Karatsuba
+1 porque es brillante y nunca se me había ocurrido que Karatsuba era un algoritmo de convolución, pero sería bueno si pudieras explicar por qué debería resolver este problema. Puedo verlo fácilmente para el caso 2x2, pero en la configuración recursiva general no veo por qué debería solucionar este problema. Me parece plausible que ni siquiera sea reparable en general, pero no lo sé.
user541686
1
@OlliNiemitalo: Bueno, la manera fácil de explicarlo es que quiero que el error relativo sea bajo en comparación con el directo O(n2)circunvolución. (Cualquier definición razonable de "bajo" funcionaría aquí ... el error relativo que obtengo con FFT es como1014que no es bajo en ninguna definición.)
user541686
1
Los dobles IEEE solo tienen una precisión de 15 a 16 dígitos decimales en el caso general. Entonces 1e-14 es un error de tamaño razonable para una secuencia de cierto número de operaciones aritméticas (a menos que elija algunos valores mágicos).
hotpaw2
1
Si alguna vez ha diseñado un sumador de coma flotante, sabrá que el exponente está determinado por el resultado de la mantisa durante la normalización. Escogiste números que producen una mantisa estrecha poco probable.
hotpaw2
3

En lugar de descartar el algoritmo de convolución rápida, ¿por qué no usar una FFT con un rango dinámico más alto?

Una respuesta a esta pregunta muestra cómo usar la biblioteca Eigen FFT con multiprecisión de refuerzo.

Mark Borgerding
fuente
2

Creo que la precisión del algoritmo Cordic puede extenderse tanto como lo desee, si es así, use un DFT entero y una longitud de palabra apropiada para su problema.

Lo mismo sería cierto con la convolución directa, use enteros muy largos.


fuente
1

La convolución de tiempo cuadrático para obtener un resultado DFT generalmente es menos precisa (puede incurrir en un ruido numérico de cuantificación más finito, debido a una estratificación más profunda de los pasos aritméticos) que el algoritmo FFT típico cuando se usan los mismos tipos aritméticos y unidades de operación.

Es posible que desee probar tipos de datos de mayor precisión (precisión cuádruple o aritmética bignum).

hotpaw2
fuente
Er, esto está usando los mismos tipos aritméticos y unidades de operación, ¿no? Claramente es más preciso. Creo que el tipo de ruido del que estás hablando no es el mismo del que estoy hablando. Las raíces de la unidad tienen una magnitud de 1, lo que significa que simplemente no pueden representar valores muy pequeños. Esto no parece estar totalmente relacionado con la cuestión de cómo se propaga el ruido a través del sistema.
user541686
Solo parece más preciso en su ejemplo porque eligió una longitud y valores donde el redondeo funcionó a su favor. Pruebe un rango de convoluciones mucho más largas con muchos más coeficientes distintos de cero con una distribución que contenga un amplio orden de magnitudes.
hotpaw2
Sin embargo, el problema que estoy tratando de resolver no tiene nada que ver con el redondeo. Esa es una cuestión diferente que estoy no trata de resolver. Los ejemplos originales que tuve fueron exactamente como lo que acaba de decir, y funcionaron bien con convolución directa, pero fueron destruidos por FFT.
user541686
El redondeo (u otros métodos de cuantificación) está involucrado en toda la aritmética de precisión finita. Algunos resultados computacionales cambian cuando se redondean, otros no o cambian menos.
hotpaw2
Nunca dije lo contrario. Lo que les acabo de decir es que el problema que estoy tratando de resolver no tiene nada que ver con el redondeo. Es un problema diferente. No me importa evitar el redondeo, pero sí me importa evitar este problema.
user541686