Algoritma Kahan: Cara Mendapatkan Perbedaan Produk yang Tepat

Saya baru-baru ini kembali ke analisis kesalahan titik apung untuk memperbaiki beberapa detail dalam edisi berikutnya buku Rendering Berbasis Fisik . Angka floating-point adalah bidang perhitungan yang menarik, penuh kejutan (baik dan buruk), serta trik rumit untuk menghilangkan kejutan yang tidak menyenangkan.

Dalam prosesnya, saya menemukan posting ini di StackOverflow , dari mana saya belajar tentang algoritma yang elegan untuk perhitungan yang tepat $$$a \ kali b-c \ kali d$ .

Tetapi sebelum melanjutkan dengan algoritma, Anda perlu memahami apa yang begitu licik dalam berekspresi $$$a \ kali b-c \ kali d$ ? Ambil $$$a = 33962.035$ , $$$b = -30438,8$ , $$$c = 41563,4$ dan $$$d = -24871.969$ . (Ini adalah nilai nyata yang saya dapatkan selama peluncuran pbrt .) Untuk nilai float 32-bit, kita dapatkan: $$$a \ kali b = -1.03376365 \ kali 10 ^ 9$ dan $$$c \ kali d = -1.03376352 \ kali 10 ^ 9$ . Kami melakukan pengurangan, dan kami mendapatkannya $$$-128$ . Tetapi jika Anda melakukan perhitungan dengan presisi ganda, dan pada akhirnya mengubahnya menjadi mengambang, Anda dapatkan $$$-75,1656$ . Apa yang terjadi

Masalahnya adalah bahwa nilai setiap pekerjaan dapat jauh melampaui garis bawah $$$-1 \ kali 10 ^ 9$ , di mana jarak antara nilai floating point yang diwakili sangat besar - 64. Artinya, saat pembulatan $$$a \ kali b$ dan $$$c \ kali d$ secara individual, ke float representable terdekat, mereka berubah menjadi angka yang merupakan kelipatan dari 64. Pada gilirannya, perbedaan mereka akan menjadi kelipatan dari 64, dan tidak akan ada harapan bahwa itu akan menjadi $$$-75,1656$ lebih dekat dari $$$-64$ . Dalam kasus kami, hasilnya bahkan lebih jauh karena bagaimana keduanya bekerja $$$-1 \ kali 10 ^ 9$ . Kita akan secara langsung menghadapi pengurangan bencana tua yang baik ¹ .

Ini adalah solusi yang lebih baik daripada ² :

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts() menghitung $$$a \ kali b-c \ kali d$ dengan cara yang menghindari kontraksi bencana. Teknik ini pertama kali dijelaskan oleh William Kahan yang legendaris dalam artikel On Cost of Floating-PointComputasi Tanpa Aritmatika Ekstra Tepat . Perlu dicatat bahwa karya Kahan menarik untuk dibaca secara keseluruhan, mereka memiliki banyak komentar tentang keadaan dunia saat ini poin mengambang, serta pertimbangan matematika dan teknis. Inilah salah satu temuannya:

Kita yang telah berjuang dengan perubahan aritmatika floating point dan "optimisasi" kompiler yang buruk dapat dengan bangga bangga dengan kemenangan dalam pertempuran ini. Tetapi jika kita meneruskan kelanjutan pertempuran ini untuk generasi masa depan, ini akan bertentangan dengan esensi seluruh peradaban. Pengalaman kami mengatakan bahwa bahasa pemrograman dan sistem pengembangan adalah sumber dari terlalu banyak kekacauan yang harus kita tangani. Terlalu banyak kesalahan dapat ditiadakan, serta beberapa "optimasi" menarik yang aman untuk bilangan bulat, tetapi kadang-kadang terbukti berakibat fatal untuk angka floating point.

Setelah membayar upeti untuk kecerdasannya, mari kita kembali ke DifferenceOfProducts() : dasar dari penguasaan fungsi ini adalah penggunaan penggandaan multiply-add, instruksi FMA ³ . Dari sudut pandang matematika, FMA(a,b,c) adalah $$$a \ kali b + c$ Karena itu, pada awalnya tampaknya operasi ini hanya berguna sebagai optimasi mikro: satu instruksi, bukan dua. Namun fma
Ini memiliki properti khusus - itu hanya membulatkan sekali.

Seperti biasa $$$a \ kali b + c$ pertama kali dihitung $$$a \ kali b$ , dan kemudian nilai ini, yang dalam kasus umum tidak dapat direpresentasikan dalam format floating point, dibulatkan ke float terdekat. Kemudian untuk nilai bulat ini ditambahkan $$$c$ , dan hasil ini kembali dibulatkan ke pelampung terdekat. FMA diimplementasikan sedemikian rupa sehingga pembulatan dilakukan hanya pada akhirnya - nilai menengah $$$a \ kali b$ mempertahankan akurasi yang cukup, jadi setelah menambahkannya $$$c$ hasil akhir akan paling mendekati nilai sebenarnya $$$a \ kali b + c$ nilai float.

Setelah berurusan dengan FMA, kami akan kembali ke DifferenceOfProducts() . Sekali lagi, saya akan menunjukkan dua baris pertama:

  float cd = c * d; float err = std::fma(-c, d, cd); 

Yang pertama menghitung nilai bulat $$$c \ kali d$ dan yang kedua ... kurangi $$$c \ kali d$ dari pekerjaan mereka? Jika Anda tidak tahu cara kerja FMA, maka Anda mungkin berpikir bahwa err akan selalu menjadi nol. Tetapi ketika bekerja dengan FMA, baris kedua sebenarnya mengekstraksi nilai kesalahan pembulatan dalam nilai yang dihitung $$$c \ kali d$ dan menyimpannya untuk err . Setelah itu, hasilnya sangat sederhana:

  float dop = std::fma(a, b, -cd); return dop + err; 

FMA kedua menghitung perbedaan antara karya-karya menggunakan FMA, melakukan pembulatan hanya di bagian paling akhir. Akibatnya, tahan terhadap pengurangan bencana, tetapi perlu bekerja dengan nilai bulat. $$$c \ kali d$ . return "memperbaiki" masalah ini, menambahkan kesalahan yang disorot di baris kedua. Dalam sebuah artikel oleh Jeannenrod et al. ditunjukkan bahwa hasilnya benar hingga 1,5 ulps (ukuran presisi unit), yang sangat bagus: FMA dan operasi floating point sederhana akurat hingga 0,5 ulps, sehingga algoritmenya hampir sempurna.

Kami menggunakan palu baru

Ketika Anda mulai mencari cara untuk menggunakan DifferenceOfProducts() , ternyata sangat berguna. Menghitung pembeda persamaan kuadrat? Call DifferenceOfProducts(b, b, 4 * a, c) ⁴ . Menghitung determinan matriks 2x2? Algoritma akan menyelesaikan masalah ini. Dalam versi pbrt berikutnya , saya menemukan sekitar 80 kegunaan untuk itu. Dari mereka semua, fungsi produk vektor adalah yang paling dicintai. Itu selalu menjadi sumber masalah, karena itu Anda harus mengangkat tangan dan menggunakan ganda dalam implementasi untuk menghindari pengurangan bencana:

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

Dan sekarang kita dapat terus bekerja dengan float dan menggunakan DifferenceOfProducts() .

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

Contoh licik dari awal posting sebenarnya adalah bagian dari karya vektor. Pada tahap tertentu, kode pbrt perlu menghitung produk vektor dari vektor $$$(33962.035, 41563.4, 7706.415)$ dan $$$(- 24871.969, -30438.8, -5643.727)$ . Ketika dihitung menggunakan float, kita akan mendapatkan vektor $$$(1552, -1248, -128)$ . (Aturan umum: jika dalam perhitungan titik apung di mana angka besar terlibat, Anda tidak mendapatkan nilai integer yang begitu besar, maka ini adalah tanda yang hampir pasti bahwa pengurangan bencana telah terjadi.)

Dengan presisi ganda, produk vektor adalah $$$(1556.0276, -1257.5151, -75.1656)$ . Kita melihatnya dengan pelampung $$$x$ terlihat normal $$$y$ sudah sangat buruk juga $$$z$ menjadi bencana yang telah menjadi motivasi untuk menemukan solusi. Dan hasil apa yang akan kita dapatkan dengan DifferenceOfProducts() dan nilai float? $$$(1556.0276, -1257.5153, -75.1656)$ . Nilai-nilai $$$x$ dan $$$z$ sesuai dengan presisi ganda, dan $$$y$ sedikit diimbangi - karenanya ulp ekstra berasal dari.

Bagaimana dengan kecepatan? DifferenceOfProducts() melakukan dua FMA, serta multiplikasi dan penambahan. Algoritma naif dapat diimplementasikan dengan satu FMA dan satu multiplikasi, yang, tampaknya, akan memakan waktu setengah. Namun dalam praktiknya, ternyata setelah mendapatkan nilai dari register, itu tidak masalah: dalam tolok ukur sintetis yang dipegang di laptop saya, DifferenceOfProducts() hanya 1,09 kali lebih mahal daripada algoritma naif. Operasi presisi ganda lebih lambat 2,98 kali.

Segera setelah Anda mengetahui tentang kemungkinan pengurangan bencana, semua jenis ekspresi yang tampak tidak bersalah dalam kode mulai tampak mencurigakan. DifferenceOfProducts() tampaknya menjadi obat yang baik untuk sebagian besar dari mereka. Mudah digunakan, dan kami tidak punya alasan khusus untuk tidak menggunakannya.

Catatan

1. Pengurangan katastropik tidak menjadi masalah ketika mengurangi jumlah dengan tanda yang berbeda atau ketika menambahkan nilai dengan tanda yang sama. Namun, ini bisa menjadi masalah saat menambahkan nilai dengan tanda yang berbeda. Artinya, jumlahnya harus dipertimbangkan dengan kecurigaan yang sama dengan perbedaan.
2. Sebagai latihan untuk pembaca, saya sarankan untuk menulis fungsi SumOfProducts() , yang memberikan perlindungan terhadap kontraksi bencana. Jika Anda ingin menyulitkan tugas, maka jelaskan mengapa di DifferenceOfProducts() , dop + err == dop , karena tanda-tanda a*b dan c*d berbeda.
3. Instruksi FMA telah tersedia di GPU selama lebih dari satu dekade, dan pada sebagian besar CPU setidaknya selama lima tahun. Mungkin perlu menambahkan flag compiler ke CPU untuk secara langsung membuatnya ketika menggunakan std::fma() ; dalam gcc dan dentang, -march=native karya -march=native .
4. Dalam format floating point IEEE, perkalian dengan kekuatan dua dilakukan persis, sehingga 4 * a tidak menyebabkan kesalahan pembulatan, kecuali terjadi overflow.