Algoritmo de convolução rápido e preciso (como FFT) para alta faixa dinâmica?

8

Parece que a convolução baseada em FFT sofre de resolução limitada de ponto flutuante devido à avaliação de tudo em torno das raízes da unidade, como você pode ver no erro de fator neste código 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]

Existem algoritmos de convolução rápida que não sofrem com esse problema?
Ou a convolução direta (tempo quadrático) é a única maneira de obter uma solução precisa?

(Se esses números pequenos são significativos o suficiente para não serem cortados é o meu ponto).

user541686
fonte
Observe que convolve()apenas chama fftconvolve()agora, se o tamanho da entrada for grande. Especifique method='direct'se você deseja direto.
endolith
@ endolith: Bom ponto! Eu aprendi isso recentemente, mas esqueci aqui.
user541686

Respostas:

5

Isenção de responsabilidade: Eu sei que esse tópico é mais antigo, mas se alguém estiver procurando por "convolução rápida precisa, alta faixa dinâmica" ou semelhante, este é um dos primeiros de apenas alguns resultados decentes. Quero compartilhar minhas idéias que obtive sobre esse tópico, para que possa ajudar alguém no futuro. Peço desculpas se posso usar os termos errados na minha resposta, mas tudo o que encontrei neste tópico é bastante vago e leva a confusão, mesmo neste tópico. Espero que o leitor entenda de qualquer maneira.

A convolução direta é mais precisa para usinar a precisão de cada ponto, ou seja, o erro relativo é geralmente aproximado ou próximo de 1.e-16 para precisão dupla para cada ponto do resultado. Cada ponto tem 16 dígitos corretos. Os erros de arredondamento podem ser significativos para convoluções atipicamente grandes e, estritamente falando, deve-se ter cuidado com o cancelamento e usar algo como a soma de Kahan e tipos de dados de precisão alta o suficiente, mas na prática o erro é quase sempre ideal.

O erro de uma convolução da FFT além dos erros de arredondamento é um erro "relativo global", o que significa que o erro em cada ponto depende da precisão da máquina e do valor de pico do resultado. Por exemplo, se o valor de pico do resultado for 2.e9, o erro absoluto em cada ponto será21091016=2107. Portanto, se um valor no resultado for muito pequeno, digamos10-9, o erro relativo nesse ponto pode ser enorme. A convolução da FFT é basicamente inútil se você precisar de pequenos erros relativos no final do seu resultado, por exemplo, se houver uma deterioração exponencial dos seus dados e se precisar de valores precisos no final. Curiosamente, se a convolução da FFT não se limita a esse erro, ela apresenta erros de arredondamento muito menores em comparação à convolução direta, pois você obviamente faz menos adições / multiplicações. Na verdade, é por isso que as pessoas costumam afirmar que a convolução da FFT é mais precisa, e elas estão quase certas em algum sentido, para que possam ser bastante inflexíveis.

Infelizmente, não há uma correção universal fácil para obter convoluções rápidas e precisas, mas, dependendo do seu problema, pode haver uma ... Encontrei duas:

Se você possui núcleos lisos que podem ser bem aproximados por um polinômio na cauda, ​​o método Multipolo rápido da caixa preta com interpolação Chebyshev pode ser interessante para você. Se o seu kernel é "bom", isso funciona perfeitamente: você obtém complexidade computacional linear (!) E precisão da precisão da máquina. Se isso se encaixa no seu problema, você deve usá-lo. No entanto, não é fácil de implementar.

Para alguns kernels específicos (acho que funções convexas, geralmente a partir de densidades de probabilidade), você pode usar um "deslocamento exponencial" para obter um erro ideal em alguma parte da cauda do resultado. Existe uma tese de doutorado e um github com uma implementação em python usando isso sistematicamente, e o autor o chama de convolução FFT precisa . Na maioria dos casos, isso não é super útil, pois ele regride de volta à convolução direta ou você pode usar a convolução da FFT de qualquer maneira. Embora o código faça isso automaticamente, o que é legal, é claro.

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

Eu olhei um pouco para o algoritmo Karatsuba (eu realmente fiz uma pequena implementação) e, para mim, parece que ele geralmente tem um comportamento de erro semelhante ao da convolução da FFT, ou seja, você recebe um erro em relação ao valor de pico do resultado. Devido à natureza de divisão e conquista do algoritmo, alguns valores na cauda do resultado realmente apresentam um erro melhor, mas não vejo uma maneira sistemática fácil de dizer quais ou, de qualquer forma, como usar essa observação. Que pena, no começo eu pensei que o Karatsuba poderia ser algo útil entre a convolução direta e a FFT. Mas não vejo casos de uso comuns em que o Karatsuba deve ser preferido em relação aos dois algoritmos comuns de convolução.

E para acrescentar à mudança exponencial que mencionei acima: há muitos casos em que você pode usá-la para melhorar o resultado de uma convolução, mas, novamente, não é uma correção universal. Na verdade, eu uso isso junto com a convolução da FFT para obter bons resultados (no caso geral de todas as entradas: no pior mesmo erro que a convolução normal da FFT, no melhor erro relativo em cada ponto da precisão da máquina). Mas, novamente, isso realmente funciona muito bem para kernels e dados específicos, mas para mim é o kernel e os dados ou um pouco exponencial em decadência.

oli
fonte
+1 bem-vindo e muito obrigado por postar isso! :)
user541686
1
Uau! Também aprendi algo e esse é um novo termo para algo que faço desde 1993. esse algoritmo de soma Kahan parece ser exatamente o mesmo que eu estava chamando de modelagem de ruído com zero na função de transferência de ruído para saída colocado à direita em DC ou o zero é colocado emz=1 no zavião. Randy Yates chamou de " economia de fração ", que é um nome genérico conciso para ele. Eu me pergunto quem é o senhor Kahan e quando isso é creditado.
Robert Bristow-johnson
2
A publicação original de Kahan parece ser de 1964.
oli
é a surpresa de hoje. Na verdade, quando um movimento @DanBoschen pediu um quebra-DSP, considerando a gama dinâmica de números de ponto flutuante, que na verdade era sobre esse mesmo conceito de adição de números muito pequenos para um número muito grande ...
Fat32
3

Um candidato é o algoritmo Karatsuba , executado emO(Nregistro23)O(N1.5849625)Tempo. Não é baseado em transformação. Há também algum código com uma explicação no arquivo de código-fonte do Music-DSP, que parece uma descoberta independente de um algoritmo semelhante.

Testar uma implementação Python do algoritmo Karatsuba (instalado por sudo pip install karatsuba) usando os números em sua pergunta mostra que, mesmo com números de ponto flutuante de 64 bits, o erro relativo é grande para um dos valores de saída:

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
fonte
2
Há um extra] no final do link para o algoritmo Karatsuba
+1 porque é brilhante e nunca me ocorreu que o Karatsuba fosse um algoritmo de convolução, mas seria bom se você pudesse explicar por que ele deveria resolver esse problema. Eu posso vê-lo facilmente no caso 2x2, mas na configuração recursiva geral, não vejo por que ele deve corrigir esse problema. Parece-me plausível que, em geral, não possa ser corrigido, mas eu não sei.
user541686
1
@OlliNiemitalo: Bem, a maneira mais fácil de explicar é que eu quero que o erro relativo seja baixo em comparação com o direto O(n2)convolução. (Qualquer definição razoável de "baixo" funcionaria aqui ... o erro relativo que estou recebendo com a FFT é como1014que não é baixo por qualquer definição).
user541686
1
As dobras IEEE têm apenas uma precisão de 15 a 16 dígitos decimais no caso geral. Portanto, 1e-14 é um erro de tamanho razoável para uma sequência de um certo número de operações aritméticas (a menos que você escolha alguns valores mágicos).
hotpaw2
1
Se você já projetou um somador de ponto flutuante, saberá que o expoente é determinado pelo resultado da mantissa durante a normalização. Você selecionou números que produzem uma mantissa estreita improvável.
hotpaw2
3

Em vez de descartar o algoritmo de convolução rápida, por que não usar uma FFT com maior faixa dinâmica?

Uma resposta a esta pergunta mostra como usar a biblioteca Eigen FFT com multiprecisão de impulso.

Mark Borgerding
fonte
2

Acredito que a precisão do algoritmo Cordic possa ser estendida até onde você desejar, se assim for, use um DFT inteiro e um comprimento de palavra apropriado ao seu problema.

O mesmo seria verdade com a convolução direta, use números inteiros muito longos.


fonte
1

A convolução quadrática no tempo para obter um resultado de DFT geralmente é menos precisa (pode incorrer em ruído numérico de quantização mais finita, devido a uma camada mais profunda de etapas aritméticas) do que o algoritmo FFT típico ao usar os mesmos tipos aritméticos e unidades de operação.

Convém tentar tipos de dados com maior precisão (precisão quádrupla ou aritmética bignum).

hotpaw2
fonte
Er, isso está usando os mesmos tipos aritméticos e unidades de operação, não é? Claramente, é mais preciso. Eu acho que o tipo de ruído que você está falando não é o mesmo que eu estou falando. As raízes da unidade têm uma magnitude de 1, o que significa que simplesmente não podem representar valores muito pequenos. Isso não parece totalmente relacionado à questão de como o ruído se propaga através do sistema.
user541686
Parece apenas mais preciso no seu exemplo, porque você escolheu um comprimento e valores em que o arredondamento funcionou a seu favor. Tente uma série de convoluções muito mais longas, com muito mais coeficientes diferentes de zero, com uma distribuição contendo uma ampla ordem de magnitudes.
precisa saber é o seguinte
O problema que estou tentando resolver não tem nada a ver com arredondamentos. Essa é uma questão diferente que não estou tentando resolver. Os exemplos originais que eu tinha eram exatamente como o que você acabou de dizer, e eles funcionaram muito bem com a convolução direta, mas foram destruídos pela FFT.
user541686
O arredondamento (ou outros métodos de quantização) está envolvido em toda aritmética de precisão finita. Alguns resultados computacionais mudam quando arredondados, outros não mudam ou mudam menos.
hotpaw2
Eu nunca reivindiquei o contrário. O que eu acabei de dizer foi o problema que estou tentando resolver não tem nada a ver com arredondamentos. É um problema diferente. Não me importo de evitar arredondamentos, mas me importo em evitar esse problema.
user541686