Perbandingan relatif dari angka floating point

10

Saya memiliki fungsi numerik yang f(x, y)mengembalikan angka floating point ganda yang mengimplementasikan beberapa rumus dan saya ingin memeriksa apakah itu benar terhadap ekspresi analitik untuk semua kombinasi parameter xdan ysaya tertarik. Apa cara yang tepat untuk membandingkan yang dihitung dan dihitung angka floating point analitik?

Katakanlah dua angka itu adalah adan b. Sejauh ini saya telah memastikan bahwa kedua kesalahan absolut ( abs(a-b) < eps) dan juga relatif ( abs(a-b)/max(abs(a), abs(b)) < eps) kurang dari eps. Dengan begitu ia akan menangkap ketidakakuratan numerik bahkan jika angkanya sekitar 1e-20.

Namun, hari ini saya menemukan masalah, nilai numerik adan nilai analitiknya badalah:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Jadi kesalahan absolut [50] adalah (jelas) kecil, tetapi kesalahan relatif [52] besar. Jadi saya pikir saya punya bug di program saya. Dengan men-debug, saya menyadari, bahwa angka-angka ini tidak normal . Karena itu, saya menulis rutin berikut untuk melakukan perbandingan relatif yang tepat:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Di mana tiny(1._dp)mengembalikan 2.22507385850720138E-308 di komputer saya. Sekarang semuanya berfungsi dan saya hanya mendapatkan 0 sebagai kesalahan relatif dan semuanya baik-baik saja. Secara khusus, kesalahan relatif di atas [52] salah, itu hanya disebabkan oleh kurang akuratnya angka denormal. Apakah implementasi rel_errorfungsi saya benar? Haruskah saya memeriksa yang abs(a-b)kurang dari kecil (= tidak normal), dan mengembalikan 0? Atau haruskah saya memeriksa beberapa kombinasi lain, seperti max(abs(a), abs(b))?

Saya hanya ingin tahu apa cara yang "tepat" itu.

floating-point Ondřej Čertík
sumber

Jawaban:

11

Anda dapat langsung memeriksa penolakan menggunakan isnormal()dari math.h(C99 atau lebih baru, POSIX.1 atau lebih baru). Di Fortran, jika modul ieee_arithmetictersedia, Anda dapat menggunakan ieee_is_normal(). Untuk lebih tepatnya tentang kesetaraan fuzzy, Anda harus mempertimbangkan representasi floating point dari denormals dan memutuskan apa yang Anda maksud agar hasilnya cukup baik.

Lebih penting lagi, untuk meyakini bahwa salah satu hasilnya benar, Anda harus yakin bahwa Anda tidak kehilangan terlalu banyak digit pada langkah perantara. Komputasi dengan denormals umumnya tidak dapat diandalkan dan harus dihindari dengan meminta algoritma Anda melakukan skala ulang secara internal. Untuk memastikan bahwa penskalaan internal Anda berhasil, saya sarankan untuk mengaktifkan pengecualian floating point menggunakan feenableexcept()dengan C99 atau ieee_arithmeticmodul di Fortran.

Meskipun aplikasi Anda dapat menangkap sinyal yang muncul pada pengecualian floating point, semua kernel yang saya coba atur ulang flag perangkat kerasnya sehingga fetestexcept()tidak mengembalikan hasil yang bermanfaat. Ketika dijalankan dengan -fp_trap, program-program PETSc akan (secara default) mencetak jejak stack ketika kesalahan floating point dinaikkan, tetapi tidak akan mengidentifikasi garis yang menyinggung. Jika Anda menjalankan debugger, debugger mempertahankan flag perangkat keras dan merusak ekspresi yang menyinggung. Anda dapat memeriksa alasan yang tepat dengan menelepon fetestexceptdari debugger di mana hasilnya adalah bitwise ATAU dari flag berikut (nilai dapat bervariasi berdasarkan mesin, lihat fenv.h; nilai-nilai ini untuk x86-64 dengan glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Jed Brown
sumber
Terima kasih atas jawabannya. Ekspresi analitik yang saya bandingkan dalam rezim asimptotik adalah exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2untuk m = 234, t = 2000. Ini menjadi nol dengan cepat saat saya bertambah m. Semua saya ingin memastikan bahwa rutin numerik saya mengembalikan angka "benar" (untuk mengembalikan nol juga baik-baik saja) ke setidaknya 12 digit signifikan. Jadi, jika perhitungan mengembalikan angka yang tidak normal, maka itu hanya nol, dan seharusnya tidak ada masalah. Jadi, hanya perbandingan rutin harus kuat terhadap ini.
Ondřej Čertík
5

Donald Knuth memiliki proposal untuk algoritma perbandingan titik mengambang dalam volume 2 "Algoritma seminumerik" dari "The Art of Computer Programming". Itu diimplementasikan dalam C oleh Th. Belding (lihat paket fcmp ) dan tersedia di GSL .

GertVdE
sumber
2
Berikut ini adalah implementasi Fortran saya: gist.github.com/3776847 , perhatikan bahwa saya harus secara eksplisit menangani angka-angka tidak normal di dalamnya. Kalau tidak, saya pikir itu cukup banyak setara dengan kesalahan relatif, satu-satunya perbedaan adalah bahwa alih-alih melakukan abs(a-b)/max(a, b) < eps, kita lakukan abs(a-b)/2**exponent(max(a, b)) < eps, yang cukup banyak hanya menjatuhkan mantissa di max(a, b), jadi menurut saya perbedaannya dapat diabaikan.
Ondřej Čertík
5

Angka denormalized yang dibulatkan secara optimal mungkin memang memiliki kesalahan relatif tinggi. (Menyiram ke nol sambil tetap menyebutnya kesalahan relatif menyesatkan.)

Tetapi mendekati nol, menghitung erros relatif tidak ada artinya.

Oleh karena itu, bahkan sebelum Anda mencapai angka denormalized, Anda mungkin harus beralih ke akurasi absolut (yaitu yang ingin Anda jamin dalam kasus ini).

yx|y−x|≤absacc+relacc∗max(|x|,|y|)

Kemudian pengguna kode Anda tahu persis seberapa akurat yang mereka miliki.

Arnold Neumaier
sumber
Apakah Anda yakin tidak ada artinya menghitung kesalahan relatif mendekati nol? Saya pikir itu tidak ada artinya hanya jika ada kehilangan keakuratan (untuk alasan apa pun). Jika misalnya ada kehilangan akurasi untuk x <1e-150 karena beberapa masalah numerik (seperti mengurangi dua angka besar), maka Anda benar. Namun dalam kasus saya angka-angkanya tampaknya akurat sampai nol, kecuali ketika angka itu mengenai angka-angka yang tidak normal. Jadi dalam kasus saya, absacc = 1e-320 atau lebih dan saya bisa memeriksa abs(a-b) < tiny(1._dp)seperti yang saya lakukan di atas.
Ondřej Čertík
@ OndřejČertík: Dalam hal ini, ganti 1e-150 dengan 1e-300 atau apa pun yang dapat Anda verifikasi. Dalam kasus apa pun yang mendekati nol Anda membuat kesalahan absolut, dan klaim kesalahan Anda harus mencerminkan hal ini daripada menyatakan kesalahan relatif menjadi nol.
Arnold Neumaier
Saya melihat. Saya dapat memverifikasi, bahwa semua berfungsi untuk angka yang lebih tinggi dari tiny(1._dp)=2.22507385850720138E-308(saya membuat kesalahan dalam komentar saya sebelumnya, ini 2e-308, bukan 1e-320). Jadi ini adalah kesalahan mutlak saya. Maka saya perlu membandingkan kesalahan relatif. Saya mengerti maksud Anda, saya pikir Anda benar. Terima kasih!
Ondřej Čertík
1
@ OndřejČertík: Untuk menemukan kesalahan relatif tambahan yang diberikan absacc, pantau maksimum . |y−x|−absaccmax(|x|,|y|)
Arnold Neumaier