Convolución, transformada rápida de Fourier y polinomios (2022)

 (alvarorevuelta.com)
1 puntos por GN⁺ 2024-07-02 | 1 comentarios | Compartir por WhatsApp
  • Si se expanden polinomios de alto grado con el método escolar, hay que multiplicar todos los pares de términos, por lo que el costo O(n²) se convierte rápidamente en un cuello de botella
  • La multiplicación de vectores de coeficientes de polinomios es igual a la convolución de señales discretas, y el resultado de [2, 3, 4] y [5, 6, 7] es [10, 27, 52, 45, 28]
  • La DFT lleva una señal discreta al dominio de la frecuencia, y la FFT calcula la misma transformación en O(n log n), marcando la diferencia con entradas grandes
  • La convolución en el dominio del tiempo se convierte en una multiplicación elemento a elemento en el dominio de la frecuencia, por lo que se puede multiplicar polinomios más rápido transformando con FFT, multiplicando y luego volviendo con IFFT
  • En grados pequeños, el costo de ida y vuelta FFT/IFFT puede anular la ventaja, pero a medida que el grado crece, el método con FFT se vuelve más eficiente

Por qué se vuelve lenta la multiplicación de polinomios

  • Un polinomio P(x) se expresa como una suma de términos formados por coeficientes a_k y potencias de la variable x
    • Por ejemplo, P(x)=5x²+2x+9 es un polinomio de grado 2
    • Según la convención de notación, el vector de coeficientes puede representarse como [5, 2, 9] o como [9, 2, 5]
  • La suma y la resta son relativamente simples, porque basta con sumar o restar los términos del mismo grado
    • En Python se puede recorrer cada coeficiente con zip(p, q) y calcular a + b o a - b
    • Si los grados son distintos, se puede usar zip_longest
  • La multiplicación requiere multiplicar cada término con los demás y luego volver a combinar los términos del mismo grado, lo que aumenta el volumen de cálculo
    • El resultado de (2x²+3x+4) × (5x²+6x+7) es 10x⁴+27x³+52x²+45x+28
    • La complejidad de este método es O(n²), y a medida que crece el grado aumenta la cantidad de multiplicaciones necesarias

Vectores de coeficientes y convolución

  • En el dominio discreto, la convolución de dos señales p y q se define como y[n]=Σ p[k]·q[n-k]
  • El cálculo consiste en invertir q, desplazarla de izquierda a derecha sobre p y sumar los productos de los elementos que se superponen
  • Las señales de ejemplo son las siguientes
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • Al invertir y desplazar q, cada coeficiente de salida se genera en el siguiente orden
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • El resultado de la convolución es y = [10, 27, 52, 45, 28]
    • Esto coincide con los coeficientes de 10x⁴+27x³+52x²+45x+28, obtenidos mediante la multiplicación de polinomios
    • Por lo tanto, la multiplicación de polinomios puede verse como una convolución de vectores de coeficientes

Transformada de Fourier y FFT

  • La transformada de Fourier convierte una señal del dominio del tiempo al dominio de la frecuencia
    • Desde la perspectiva del tiempo, la señal se ve como valores en instantes específicos
    • Desde la perspectiva de la frecuencia, la señal se interpreta como la suma de distintas frecuencias de oscilación
  • Las frecuencias de oscilación se representan con senos y cosenos, cada uno con coeficiente y fase
  • Si se aplica FFT a una onda sinusoidal pura de 5 Hz, en el dominio de la frecuencia aparece como un delta en la posición de 5 Hz
    • Esto muestra que la onda sinusoidal en el dominio del tiempo puede representarse como un único seno de 5 Hz
  • Los términos relacionados se distinguen así
    • Fourier Transform(FT): transformada de Fourier definida en el dominio continuo
    • Discrete Fourier Transform(DFT): transformada de Fourier definida para señales discretas
    • Fast Fourier Transform(FFT): algoritmo que calcula la DFT en O(n log n) en lugar de O(n²)
  • La DFT convierte una señal de tiempo discreto x[n] en X[k] en el dominio de la frecuencia
    • Cada X[k] se calcula multiplicando y sumando las muestras de entrada por números complejos que representan una frecuencia específica

Convertir en multiplicación en el dominio de la frecuencia

  • La ventaja central de la DFT y del dominio de la frecuencia es que la convolución puede convertirse en multiplicación elemento a elemento
    • Convolucionar dos señales en el dominio del tiempo equivale a multiplicarlas en el dominio de la frecuencia
    • La multiplicación se puede calcular más rápido que la convolución
  • El procedimiento para multiplicar polinomios rápidamente es el siguiente
    • Transformar los polinomios al dominio de la frecuencia con FFT: O(n log n)
    • Multiplicar elemento a elemento en el dominio de la frecuencia: O(n)
    • Transformar el resultado de vuelta al dominio del tiempo con IFFT: O(n log n)
  • En conjunto, al usar FFT se puede realizar la multiplicación de polinomios con complejidad O(n log n)
  • En polinomios grandes, es más rápido que la multiplicación escolar O(n²)

Implementación en Python y benchmark

  • multiply_naive multiplica todos los pares de coeficientes con un doble bucle y los suma en la posición de resultado i + j
    • La longitud del resultado es len(p) + len(q) - 1
    • Su complejidad es O(n²)
  • multiply_fft realiza la multiplicación de coeficientes con base en FFT/IFFT
    • Calcula una longitud que sea potencia de 2 y al menos len(p) + len(q) - 1, para que quepa el resultado
    • Rellena ambas entradas con np.pad
    • Multiplica elemento a elemento los valores transformados con np.fft.fft
    • Vuelve con np.fft.ifft y luego redondea la parte real para convertirla en coeficientes enteros
  • Con la entrada de ejemplo p = [2, 3, 4], q = [5, 6, 7], ambos métodos devuelven [10, 27, 52, 45, 28]
  • En el benchmark se compara el método con FFT con multiply_convolve, que usa np.convolve en lugar de multiply_naive
    • Esto se debe a que multiply_naive es lento por los bucles de Python, por lo que es difícil compararlo directamente con el método FFT que usa np.fft.fft
    • np.convolve realiza la misma operación en código C de bajo nivel
  • El grado se incrementa en el rango range(1, 30000, 1000), y para cada grado se generan dos polinomios con coeficientes aleatorios entre 1 y 999999
    • Para cada método se mide el tiempo promedio con n_runs = 5
    • En grados bajos, el método FFT puede no ser ventajoso por el costo de las transformaciones de ida y vuelta FFT/IFFT
    • A medida que aumenta el grado, el método FFT muestra resultados mucho más eficientes

Lecturas relacionadas

1 comentarios

 
GN⁺ 2024-07-02
Comentarios de Hacker News

  • Lo que siempre me molesta de estas explicaciones es que normalmente se olvidan del error numérico
    No se puede abstraer la multiplicación de coeficientes simplemente como si fuera de “tiempo constante”. Si vas a hacer eso, entonces también podrías abstraer toda la multiplicación desde el principio. Si se toma en cuenta la precisión numérica, está más cerca de O(n (log n)^3) [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • El límite de error citado en ese artículo es excesivamente pesimista. La edición más reciente de Knuth ya incluye el límite correcto, porque yo se lo hice saber
    • Ojalá se pudiera usar la aritmética basada en cuaterniones que aparece en la publicación del OP para reducir o incluso eliminar el error de multiplicación [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 los coeficientes son enteros, se puede obtener un resultado exacto con una NTT usando un módulo suficientemente grande, y sobre todo en hardware el tiempo de multiplicación también puede ser más rápido
    • Por eso se distingue entre ciencias de la computación e ingeniería de software :)

  • Con este método se pueden multiplicar números largos entre sí. La idea clave es que la multiplicación de polinomios es equivalente a la multiplicación usual de números largos sin hacer acarreos
    Por ejemplo, si tienes un número de 1000 dígitos, tomas cada dígito como coeficiente de un polinomio con 1000 elementos. Luego puedes multiplicar esos polinomios con el método FFT explicado en el artículo. Para convertir el resultado de nuevo en un número, hay que procesar los acarreos. Si algún elemento es mayor que 10, se pasa el exceso al siguiente dígito, y luego se convierten los coeficientes en un número
    Esa es la idea básica, aunque hay sutilezas en la precisión necesaria para los acarreos y en garantizar que redondear el resultado de la FFT al entero más cercano siga siendo correcto. Así es como GMP, la biblioteca más representativa de este campo, multiplica números grandes

    • Como dijiste, tiene sentido porque un número decimal puede representarse como un polinomio evaluado en x=10. Por ejemplo, 983 = 9x^2 + 8x + 3, o sea [9, 8, 3]
      Me da curiosidad qué tan grande tiene que ser un número para que esto realmente valga la pena, y en qué se usa

  • Si todavía no lo has visto, vale la pena ver este video
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    Es realmente excelente para derivar el algoritmo FFT a partir de la multiplicación de polinomios. Lo vuelvo a ver como cada 6 meses

  • La propiedad de la FFT de que “la convolución es multiplicación punto a punto” también se cumple en cualquier grupo multiplicativo cíclico. Para una derivación más algebraica, ver https://www.sciencedirect.com/science/article/pii/S002200007...
    A veces a esto se le llama “FFT armónica”, y también hay FFT no armónicas: la [LCH14] “additive NTT” sobre GF(2^n), la circle FFT de [HLP24] sobre la variedad unitaria X^2+Y^2=1 en campos finitos, y ecfft de [BCKL21] sobre torres de isogenias de curvas elípticas
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • ¿Quién habrá sido la primera persona en proponer usar la FFT para una multiplicación de polinomios más rápida?
    Me dio curiosidad hace poco y me puse a investigar; aunque no soy muy bueno siguiendo citas, logré remontarme hasta el artículo de 1995 de David Eppstein [0]. Ahí se usa para resolver eficientemente el problema de suma de subconjuntos tras actualizaciones incrementales. Seguro que en el TAOCP de Knuth debe haber algo incluso anterior
    También me impactó bastante saber que la suma de subconjuntos exacta permitiendo repeticiones puede resolverse en tiempo subexponencial con multiplicación de polinomios por FFT [1]. Lo importante es que este algoritmo es O(N log N), donde N no es el tamaño del conjunto sino el valor máximo del elemento, así que no es algo como un contraejemplo a P ≠ NP
    [0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
    [1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...

  • Creo que todo el aprendizaje automático consiste en resolver ecuaciones de convolución
    Este artículo lo trata en el contexto del aprendizaje por refuerzo https://arxiv.org/abs/1712.06115, pero la mayoría de los enfoques encajan dentro de ese paradigma

    • ¿No significa eso, básicamente, métodos de kernel?

  • Justo implementé un algoritmo (matrix profile) que usa FFT para calcular productos internos sobre un gran conjunto de subsecuencias de series temporales. La longitud n de la serie puede llegar a cientos de millones
    Con la convolución rápida vía FFT, el tiempo de cómputo baja de O(n) a O(log n), y a esta escala la mejora de velocidad es enorme. Si además usas GPU, se acelera todavía más: algo como procesar 10 millones de puntos de datos en 0.1 segundos desde una laptop

  • El “truco” clave de esta operación parece ser esta observación:

    En otras palabras, hacer la convolución de dos señales en el dominio del tiempo es lo mismo que multiplicar esas dos señales en el dominio de la frecuencia.
    Es un muy buen texto que descompone una idea compleja en pasos mucho más pequeños, de modo que incluso yo, que no soy fuerte en matemáticas, pude entenderlo más o menos. Pero, ¿me habré perdido un paso intermedio? ¿O lo dejaron como ejercicio para que el lector lo investigue? Para ese punto yo ya estaba usando al máximo toda mi capacidad matemática, y me dio un poco la sensación de “y ahora dibuja el maldito resto del búho”. ¿Soy el único? Igual, el texto me pareció muy bueno

    • No sé si ayude, pero: la multiplicación de dos polinomios como se enseña en la escuela en realidad es una convolución
      Existe la propiedad de que “hacer la convolución de dos señales en el dominio del tiempo es lo mismo que multiplicar esas dos señales en el dominio de la frecuencia”, y la FFT permite transformar del dominio del tiempo al dominio de la frecuencia. Entonces pasas los polinomios al dominio de la frecuencia con FFT, y ahí solo multiplicas. Eso es más rápido que hacer la convolución. Me pregunto si eso aclara el paso faltante; si todavía falta algo, puedo actualizar el texto

  • Entonces, ¿la factorización de enteros es una deconvolución discreta? Me pregunto si, al poner lado a lado la representación por FFT, o sea la multiplicación punto a punto, y el tableau, es decir la multiplicación larga normal con suma de acarreos, la simetría se rompe y se puede obtener suficiente información para un algoritmo rápido

  • Claro, la multiplicación ingenua de polinomios es lenta respecto al grado del polinomio. Pero, ¿cuándo realmente hace falta trabajar con dos polinomios de grado 100?
    Por eso da la impresión de que en los sistemas de álgebra computacional no se usan estos métodos.

    • Los sistemas de álgebra computacional, por ejemplo chebfun de Matlab, convierten funciones arbitrarias en polinomios de grado mayor que 100 para encontrar raíces, óptimos, etc. con más facilidad.
    • En corrección de errores y procesamiento de señales es muy común.
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... es un ejemplo.
    • Quise hacer ingeniería inversa de los parámetros del checksum CRC de archivos grandes, así que hice un programa[1] que convierte archivos en polinomios GF(2) de millones de grados y calcula el máximo común divisor. Sin multiplicación basada en FFT, sería imposible hacerlo en un tiempo razonable.
      [1]: https://github.com/8051enthusiast/delsum
    • Esta perspectiva de convolución y los kernels rápidos de GPU para FFT se usaron en algunos modelos de espacio de estados anteriores a Mamba para modelar secuencias largas; aquí el polinomio es la secuencia de entrada.
      En las entradas del blog de Hazy Research de 2020 a 2023 hay mucha información sobre este enfoque.
    • Ver https://news.ycombinator.com/item?id=40306339
      "(...) en investigación de física he trabajado con expresiones de casi 1 terabyte de longitud y con más de 100 millones de términos"