Algoritma konvolusi cepat & akurat (seperti FFT) untuk rentang dinamis tinggi?

Tampaknya konvolusi berbasis FFT menderita resolusi floating-point terbatas karena mengevaluasi segala sesuatu di sekitar akar persatuan, seperti yang dapat Anda lihat di $10^{14}$-faktor kesalahan dalam kode Python ini:

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]

Adakah algoritma konvolusi cepat yang tidak mengalami masalah ini?
Atau apakah konvolusi langsung (waktu kuadratik) satu-satunya cara untuk mendapatkan solusi yang akurat?

(Apakah angka sekecil itu cukup signifikan untuk tidak terpotong adalah poin saya.)

Penafian: Saya tahu topik ini lebih lama, tetapi jika seseorang mencari "konvolusi cepat rentang dinamis tinggi" atau serupa ini adalah salah satu yang pertama dari hanya beberapa hasil yang layak. Saya ingin berbagi wawasan saya tentang topik ini sehingga mungkin membantu seseorang di masa depan. Saya minta maaf jika saya mungkin menggunakan istilah yang salah dalam jawaban saya, tetapi semua yang saya temukan pada topik ini agak kabur dan menyebabkan kebingungan bahkan di utas ini. Saya harap pembaca akan mengerti juga.

Konvolusi langsung sebagian besar akurat ke presisi mesin untuk setiap titik, yaitu kesalahan relatif biasanya kira-kira atau mendekati 1.e-16 untuk presisi ganda untuk setiap titik hasil. Setiap titik memiliki 16 digit yang benar. Namun kesalahan pembulatan dapat signifikan untuk konvolusi besar yang tidak tipikal, dan secara tegas seseorang harus berhati-hati dengan pembatalan dan menggunakan sesuatu seperti penjumlahan Kahan dan tipe data presisi yang cukup tinggi, tetapi dalam praktiknya kesalahan tersebut hampir selalu optimal.

Kesalahan konvolusi FFT selain dari kesalahan pembulatan adalah kesalahan "relatif global", yang berarti kesalahan di setiap titik tergantung pada presisi mesin dan nilai puncak hasilnya. Misalnya jika nilai puncak dari hasilnya adalah 2.e9, maka kesalahan absolut di setiap titik adalah$2\cdot10^9\cdot10^{-16} = 2\cdot10^{-7}$. Jadi jika nilai dalam hasil seharusnya sangat kecil, katakanlah$10^{-9}$, kesalahan relatif pada titik itu bisa sangat besar. Konvolusi FFT pada dasarnya tidak berguna jika Anda memerlukan kesalahan relatif kecil di bagian hasil Anda, misalnya Anda memiliki peluruhan data yang agak eksponensial dan membutuhkan nilai yang akurat di bagian ekor. Menariknya jika konvolusi FFT tidak dibatasi oleh kesalahan itu, ia memiliki kesalahan pembulatan yang jauh lebih kecil dibandingkan dengan konvolusi langsung, karena Anda jelas melakukan penambahan / perkalian yang lebih sedikit. Inilah sebenarnya mengapa orang sering mengklaim konvolusi FFT lebih akurat, dan mereka dalam beberapa hal hampir benar, sehingga mereka bisa sangat bersikeras.

Sayangnya tidak ada perbaikan universal yang mudah untuk mendapatkan konvolusi yang cepat dan akurat, tetapi tergantung pada masalah Anda mungkin ada satu ... Saya telah menemukan dua:

Jika Anda memiliki kernel halus yang dapat didekati dengan baik oleh polinomial di bagian ekor, maka Metode Multipole Cepat kotak hitam dengan interpolasi Chebyshev mungkin menarik bagi Anda. Jika kernel Anda "baik" ini berfungsi dengan sangat sempurna: Anda mendapatkan kompleksitas komputasi linier (!) Dan akurasi presisi mesin. Jika ini cocok dengan masalah Anda, Anda harus menggunakannya. Namun tidak mudah untuk diimplementasikan.

Untuk beberapa kernel tertentu (saya pikir fungsi cembung, biasanya dari kepadatan probabilitas) Anda dapat menggunakan "pergeseran eksponensial" untuk mendapatkan kesalahan optimal di beberapa bagian ekor hasil. Ada tesis PHD dan github dengan implementasi python menggunakan itu secara sistematis, dan penulis menyebutnya konvolusi FFT akurat . Namun dalam kebanyakan kasus ini tidak terlalu berguna, karena ia akan kembali ke konvolusi langsung atau Anda tetap dapat menggunakan konvolusi FFT. Meskipun kode melakukannya secara otomatis, tentu saja itu bagus.

-------------------- EDIT: --------------------

Saya melihat sedikit pada algoritma Karatsuba (saya sebenarnya membuat implementasi kecil), dan bagi saya sepertinya memiliki perilaku kesalahan yang sama seperti konvolusi FFT, yaitu Anda mendapatkan kesalahan relatif terhadap nilai puncak hasil. Karena sifat membagi dan menaklukkan algoritma beberapa nilai di ekor hasil sebenarnya memiliki kesalahan yang lebih baik, tapi saya tidak melihat cara sistematis yang mudah untuk mengetahui mana atau dalam hal apa pun bagaimana menggunakan pengamatan ini. Sayang sekali, pada awalnya saya pikir Karatsuba mungkin sesuatu yang berguna di antara konvolusi langsung dan FFT. Tapi saya tidak melihat kasus penggunaan umum di mana Karatsuba harus lebih disukai daripada dua algoritma konvolusi umum.

Dan untuk menambah pergeseran eksponensial yang saya sebutkan di atas: Ada banyak kasus di mana Anda dapat menggunakannya untuk meningkatkan hasil konvolusi, tetapi sekali lagi ini bukan perbaikan universal. Saya benar-benar menggunakan ini bersama dengan konvolusi FFT untuk mendapatkan hasil yang cukup bagus (dalam kasus umum untuk semua input: pada kesalahan terburuk yang sama seperti konvolusi FFT normal, pada kesalahan relatif terbaik di setiap titik ke presisi mesin). Tetapi sekali lagi, ini hanya benar-benar berfungsi dengan baik untuk kernel dan data tertentu, tetapi bagi saya kernel dan data atau agak eksponensial dalam peluruhan.

Salah satu kandidat adalah algoritma Karatsuba , yang berjalan di$O\big(N^{\log_23}\big) \approx O\big(N^{1.5849625}\big)$waktu. Itu tidak berbasis transformasi. Ada juga beberapa kode dengan penjelasan di Music-DSP Source Code Archive, yang terlihat seperti penemuan independen dari algoritma yang sama.

Menguji implementasi Python dari algoritma Karatsuba (diinstal oleh sudo pip install karatsuba) menggunakan angka-angka dalam pertanyaan Anda menunjukkan bahwa bahkan dengan angka floating point 64-bit kesalahan relatif besar untuk salah satu nilai output:

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

yang mencetak:

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

Daripada membuang algoritma konvolusi cepat, mengapa tidak menggunakan FFT dengan rentang dinamis yang lebih tinggi?

Jawaban untuk pertanyaan ini menunjukkan bagaimana menggunakan perpustakaan Eigen FFT dengan meningkatkan multiprecision.

Saya percaya bahwa ketepatan algoritma Cordic dapat diperpanjang sejauh yang Anda inginkan, jika demikian gunakan DFT integer dan panjang kata yang sesuai dengan masalah Anda.

Hal yang sama berlaku untuk konvolusi langsung, gunakan bilangan bulat yang sangat panjang.

Konvolusi waktu kuadratik untuk mendapatkan hasil DFT biasanya kurang akurat (dapat menimbulkan derau numerik kuantisasi yang lebih terbatas, karena pelapisan langkah aritmatika yang lebih dalam) daripada algoritme FFT tipikal ketika menggunakan tipe aritmatika dan unit operasi yang sama.

Anda mungkin ingin mencoba tipe data dengan presisi lebih tinggi (presisi quad atau bignum aritmatika).