Perhitungan yang kuat dari rata-rata dua angka dalam floating-point?

15

Biarkan x, ymenjadi dua angka floating-point. Apa cara yang benar untuk menghitung artinya?

Cara naif (x+y)/2dapat menghasilkan luapan saat xdan yterlalu besar. Saya pikir 0.5 * x + 0.5 * ymungkin lebih baik, tetapi melibatkan dua perkalian (yang mungkin tidak efisien), dan saya tidak yakin apakah itu cukup baik. Apakah ada cara yang lebih baik?

Gagasan lain yang saya mainkan adalah (y/2)(1 + x/y)jika x<=y. Tetapi sekali lagi, saya tidak yakin bagaimana menganalisis ini dan membuktikan bahwa itu memenuhi persyaratan saya.

Selain itu, saya perlu jaminan bahwa rata-rata yang dihitung akan >= min(x,y)dan <= max(x,y). Seperti yang ditunjukkan dalam jawaban Don Hatch , mungkin cara yang lebih baik untuk mengajukan pertanyaan ini adalah: Apa implementasi dari rata-rata dua angka yang selalu memberikan hasil akurat yang paling mungkin? Yaitu, jika xdan yapakah angka floating-point, bagaimana cara menghitung angka floating-point terdekat (x+y)/2? Dalam hal ini, rata-rata yang dihitung secara otomatis >= min(x,y)dan <= max(x,y). Lihat jawaban Don Hatch untuk detailnya.

Catatan: Prioritas saya adalah akurasi yang kuat. Efisiensi dapat dihabiskan. Namun, jika ada banyak algoritma yang kuat dan akurat, saya akan memilih yang paling efisien.

becko
sumber
(+1) Pertanyaan yang menarik, tidak mengejutkan sepele.
Kirill
1
Di masa lalu, nilai floating point dihitung dan disimpan dalam bentuk presisi yang lebih tinggi untuk hasil antara. Jika a + b (64-bit dobel) menghasilkan hasil menengah 80 bit dan ini adalah apa yang dibagi 2, Anda tidak perlu khawatir meluap. Kehilangan presisi kurang jelas.
JDługosz
Solusi untuk ini tampaknya relatif sederhana ( saya menambahkan jawaban ). Masalahnya adalah saya seorang programmer dan bukan ahli ilmu komputer, jadi apa yang saya lewatkan yang membuat pertanyaan ini jauh lebih sulit?
IQAndreas
Jangan khawatir tentang biaya perkalian dan pembagian oleh dua; kompiler Anda akan mengoptimalkannya untuk Anda.
Federico Poloni

Jawaban:

18

Saya pikir Akurasi dan Stabilitas Algoritma Numerik Higham membahas bagaimana orang dapat menganalisis jenis masalah ini. Lihat Bab 2, terutama latihan 2.8.

Dalam jawaban ini saya ingin menunjukkan sesuatu yang tidak benar-benar dibahas dalam buku Higham (sepertinya tidak terlalu dikenal luas, dalam hal ini). Jika Anda tertarik untuk membuktikan properti dari algoritma numerik sederhana seperti ini, Anda dapat menggunakan kekuatan pemecah SMT modern ( Satisfiability Modulo Theories ), seperti z3 , menggunakan paket seperti sbv di Haskell. Ini agak lebih mudah daripada menggunakan pensil dan kertas.

Misalkan saya diberi nilai , dan saya ingin tahu apakah z = ( x + y ) / 2 memenuhi x z y . Berikut kode Haskell0xyz=(x+y)/2xzy

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

akan membiarkan saya melakukan ini secara otomatis . Berikut test1 funadalah proposisi bahwa untuk semua pelampung terbatas x , y dengan 0 x y .xfkamun(x,y)yx,y0xy

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

Itu meluap. Misalkan saya sekarang mengambil formula Anda yang lain: z=x/2+y/2

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

Tidak berfungsi (karena underflow bertahap: , yang mungkin tidak intuitif karena semua aritmatika menjadi basis-2).(x/2)×2x

Sekarang coba :z=x+(y-x)/2

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

Bekerja! Ini Q.E.D.adalah bukti bahwa test1properti berlaku untuk semua pelampung seperti yang didefinisikan di atas.

Bagaimana dengan yang sama, tetapi terbatas pada (bukan 0 x y )?xy0xy

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

Oke, jadi jika meluap, bagaimana dengan z = x + ( y / 2 - x / 2 ) ?y-xz=x+(y/2-x/2)

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

Jadi sepertinya di antara rumus yang saya coba di sini, tampaknya berfungsi (dengan bukti juga). Pendekatan SMT solver bagi saya merupakan cara yang jauh lebih cepat untuk menjawab kecurigaan tentang rumus titik-mengambang sederhana daripada melalui analisis kesalahan titik-mengambang dengan pensil dan kertas.x+(y/2-x/2)

Akhirnya, tujuan akurasi dan stabilitas seringkali bertentangan dengan tujuan kinerja. Untuk kinerja, saya tidak benar-benar melihat bagaimana Anda dapat melakukan lebih baik daripada , terutama karena kompiler masih akan melakukan tugas berat menerjemahkan ini menjadi instruksi mesin untuk Anda.(x+y)/2

xx+(y/2-x/2)ySFloatSDouble

-ffast-math(x+y)/2

PPPS Saya terbawa sedikit melihat hanya pada ekspresi aljabar sederhana tanpa persyaratan. Don Hatch 's rumus secara ketat baik.

Kirill
sumber
2
Tahan; apakah Anda mengklaim bahwa jika x <= y (terlepas dari apakah x> = 0 atau tidak) maka x + (y / 2-x / 2) adalah cara yang baik untuk melakukannya? Menurut saya itu tidak mungkin benar, karena itu memberikan jawaban yang salah dalam kasus berikut ketika jawaban itu benar-benar mewakili: x = -1, y = 1 + 2 ^ -52 (angka keterwakilan terkecil yang lebih besar dari 1), dalam hal ini jawabannya adalah 2 ^ -53. Konfirmasi dengan python: >>> x = -1.; y = 1.+2.**-52; print `2**-53`, `(x+y)/2.`, `x+(y/2.-x/2.)`
Don Hatch
2
x(x+y)/2yx,y(x+y)/2akan selalu menghasilkan hasil yang dibulatkan dengan benar dengan kesalahan relatif kecil. Saya sendiri lebih suka(x+y)/2.
Kirill
8

Pertama, perhatikan bahwa jika Anda memiliki metode yang memberikan jawaban paling akurat dalam semua kasus, maka itu akan memenuhi kondisi yang Anda butuhkan. (Perhatikan bahwa saya mengatakan sebuah jawaban yang paling akurat daripada yang jawabannya paling akurat, karena mungkin ada dua pemenang.) Bukti: Jika, sebaliknya, Anda punya jawaban yang akurat-sebagai-mungkin bahwa tidak tidak memenuhi kondisi yang diperlukan, yang berarti answer<min(x,y)<=max(x,y)(dalam hal min(x,y)ini jawaban yang lebih baik, kontradiksi), atau min(x,y)<=max(x,y)<answer(dalam hal inimax(x,y) ini jawaban yang lebih baik, kontradiksi).

Jadi saya pikir itu berarti pertanyaan Anda bermuara pada menemukan jawaban yang paling akurat. Dengan asumsi aritmetika IEEE754 secara keseluruhan, saya mengusulkan yang berikut:

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

Argumen saya bahwa ini memberikan jawaban yang paling akurat adalah analisis kasus yang agak membosankan. Ini dia:

  • Kasus max(abs(x),abs(y)) >= 1.:

    • Subcase baik x maupun y tidak dinormalisasi: Dalam hal ini jawaban yang dikomputasi x/2.+y/2.memanipulasi mantra yang sama dan oleh karena itu memberikan jawaban yang sama persis seperti perhitungan (x+y)/2akan menghasilkan jika kita mengasumsikan eksponen diperluas untuk mencegah luapan. Jawaban ini mungkin tergantung pada mode pembulatan tetapi dalam hal apa pun itu dijamin oleh IEEE754 untuk menjadi jawaban terbaik (dari kenyataan bahwa yang dihitung x+ydijamin menjadi perkiraan terbaik untuk matematika x + y, dan pembagian dengan 2 tepat dalam hal ini kasus).
    • Huruf x didenormalkan (dan sebagainya abs(y)>=1):

      answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

    • Subkotak y didenormalkan (dan sebagainya abs(x)>=1): analog.

  • Kasus max(abs(x),abs(y)) < 1.:
    • Subcase yang dikomputasi x+yadalah non-denormalized atau denormalized-and- "even": Meskipun dikomputasi x+ymungkin tidak tepat, itu dijamin oleh IEEE754 untuk menjadi perkiraan terbaik untuk matematika x + y. Dalam hal ini pembagian selanjutnya dengan 2 dalam ekspresi (x+y)/2.adalah tepat, sehingga jawaban yang dihitung (x+y)/2.adalah perkiraan terbaik untuk matematika (x + y) / 2.
    • Subcase yang dikomputasi didenormalisasi x+ydan "ganjil": Dalam hal ini tepat salah satu dari x, y juga harus didenormalkan-dan- "ganjil", yang berarti yang lain dari x, y didenormalkan dengan tanda yang berlawanan, sehingga yang dihitung x+yadalah persis matematis x + y, dan yang dihitung (x+y)/2.dijamin oleh IEEE754 menjadi perkiraan terbaik untuk matematika (x + y) / 2.
Don Hatch
sumber
Saya menyadari ketika saya mengatakan "didenormalisasi", saya benar-benar bermaksud sesuatu yang lain - yaitu, angka yang dekat satu sama lain dengan angka, yaitu kisaran angka yang kira-kira dua kali lebih besar dari kisaran angka yang didenormalkan, yaitu 8 kutu pertama atau lebih dalam diagram di en.wikipedia.org/wiki/Denormal_number . Intinya adalah, yang "aneh" dari ini adalah satu-satunya angka yang membaginya dengan dua tidak tepat. Saya perlu menyusun kembali bagian dari jawaban ini untuk membuatnya jelas.
Don Hatch
Anda dapat menyederhanakan analisis Anda dengan mencatat bahwa, jika tidak ada overflow / underflow, ia selalu berpendapat demikian fl(Haihal(x,y))=Haihal(x,y)(1+δ) dimana |δ|kamu, dan bahwa pembagian dengan 2 adalah tepat untuk angka non-normal. Karena itu berarti keduanyax/2+y/2 dan (x+y)/2selalu dibulatkan dengan benar, tidak ada over- / underflow, yang tersisa adalah tidak menunjukkan over- / underflow, yang mudah.
Kirill
@ Kirill aku agak tersesat ... dari mana asalmu? Juga saya tidak berpikir itu benar bahwa "pembagian dengan 2 adalah tepat untuk angka-angka non-normal" ... ini adalah hal yang sama dengan yang saya lakukan, dan tampaknya agak canggung untuk mencoba memperbaikinya. Pernyataan yang tepat adalah sesuatu yang lebih seperti "x / 2 tepat selama abs (x) setidaknya dua kali angka subnormal terbesar" ... argh, canggung!
Don Hatch
3

Untuk format titik-mengambang biner IEEE-754, dicontohkan oleh binary64perhitungan (presisi ganda), S. Boldo secara resmi membuktikan bahwa algoritma sederhana yang ditunjukkan di bawah ini memberikan rata-rata bulat yang benar.

Sylvie Boldo, "Verifikasi formal program menghitung rata-rata titik mengambang." Dalam Konferensi Internasional tentang Metode Teknik Formal , hlm. 17-32. Springer, Cham, 2015. ( konsep online )

Karena pembagian oleh dua akurat dalam aritmatika titik-mengambang biner, asalkan aliran bawah tidak terjadi , tampaknya secara intuitif jelas bahwa dengan memilih salah satu dari dua rumus(x+y)/2 dan x/2+y/2yang sesuai (berdasarkan besarnya input) sekali harus mencapai rata-rata yang dibulatkan secara akurat. Kertas Boldo menunjukkan bahwa untuk IEEE-754 binary64ada titik perpindahanC[2-967,2970]akan cukup. Orang mungkin memilihC sehingga memberikan kinerja terbaik untuk use case tertentu.

Ini menghasilkan ISO-C99kode teladan berikut :

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

Dalam pekerjaan tindak lanjut baru-baru ini, S. Boldo dan rekan penulis menunjukkan bagaimana untuk mencapai hasil terbaik untuk format floating-point desimal IEEE-754 dengan memanfaatkan operasi FMA multiply-add (FMA) dan presisi yang terkenal. blok bangunan penggandaan (TwoSum):

Sylvie Boldo, Florian Faissole, dan Vincent Tourneur, "Algoritma yang Terbukti Secara Resmi untuk Menghitung Rata-Rata yang Benar dari Angka Titik Apung Desimal." Dalam Simposium IEEE 25 tentang Aritmatika Komputer (ARITH 25) , Juni 2018, hlm. 69-75. ( konsep online )

njuffa
sumber
2

Meskipun ini mungkin bukan kinerja-efisien super-bijaksana, ada cara yang sangat sederhana untuk (1) memastikan tidak ada angka yang lebih besar dari salah satu xatau y(tidak ada luapan) dan (2) menjaga titik mengambang sebagai "akurat" seperti mungkin (dan (3) , sebagai bonus tambahan, meskipun pengurangan digunakan, tidak ada nilai yang akan disimpan sebagai angka negatif.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

Bahkan, jika Anda benar - benar ingin mendapatkan akurasi, Anda bahkan tidak perlu melakukan pembagian di tempat; cukup kembalikan nilai-nilai min(x, y)dan differenceyang dapat Anda gunakan untuk menyederhanakan secara logis atau memanipulasi nanti.

IQAndreas
sumber
Apa yang saya coba cari tahu sekarang adalah bagaimana membuat jawaban yang sama ini berfungsi dengan lebih dari dua item , sambil menjaga semua variabel tetap lebih rendah dari yang terbesar dari angka, dan menggunakan hanya satu operasi divisi untuk menjaga akurasi.
IQAndreas
@becko Yup, Anda akan melakukan pembagian setidaknya dua kali. Juga, contoh yang Anda berikan akan membuat jawaban salah. Bayangkan rata-rata 2,4,9, tidak sama dengan rata-rata 3,9.
IQAndreas
Anda benar, rekursi saya salah. Saya tidak yakin bagaimana cara memperbaikinya sekarang, tanpa kehilangan presisi.
becko
Bisakah Anda membuktikan bahwa ini memberikan hasil yang paling akurat? Yaitu, jika xdan yapakah floating-point, perhitungan Anda menghasilkan floating-point terdekat (x+y)/2?
becko
1
Tidakkah ini akan meluap ketika x, y adalah angka yang paling mudah diungkapkan dan terbesar?
Don Hatch
1

Konversikan ke precission yang lebih tinggi, tambahkan nilainya di sana dan konversi kembali.

Seharusnya tidak ada overflow dalam precisi yang lebih tinggi dan jika keduanya berada dalam kisaran floating point yang valid, jumlah yang dihitung harus di dalam juga.

Dan itu harus di antara mereka, kasus terburuk hanya setengah dari jumlah yang lebih besar jika precisi tidak mencukupi.

leeroy
sumber
Ini adalah pendekatan brute force. Mungkin berhasil, tetapi saya sedang mencari analisis yang tidak memerlukan presisi menengah yang lebih tinggi. Juga, dapatkah Anda memperkirakan berapa banyak presisi menengah yang lebih tinggi diperlukan? Bagaimanapun, jangan hapus jawaban ini (+1), saya tidak akan menerimanya sebagai jawabannya.
becko
1

Secara teoritis, x/2dapat dihitung dengan mengurangi 1 dari mantissa.

Namun, sebenarnya menerapkan operasi bitwise seperti ini tidak selalu langsung, terutama jika Anda tidak tahu format angka floating point Anda.

Jika Anda dapat melakukan ini, seluruh operasi dikurangi menjadi 3 tambahan / pengurangan, yang seharusnya merupakan peningkatan yang signifikan.

Roland Heath
sumber
0

Saya berpikir sepanjang jalan yang sama dengan @Rand Heath tetapi belum bisa berkomentar, inilah pendapat saya:

x/2dapat dihitung dengan mengurangi 1 dari eksponen (bukan mantissa, mengurangi 1 dari mantissa mengurangi 2^(value_of_exponent-length_of_mantissa)dari nilai keseluruhan).

Tanpa batasan kasus umum, mari kita asumsikan x < y. (Jika x > y, beri label ulang variabel. Jika x = y, (x+y) / 2sepele.)

  • Berubah (x+y) / 2menjadi x/2 + y/2, yang dapat dilakukan oleh dua pengurangan integer (oleh satu dari eksponen)
    • Namun ada batas bawah pada eksponen tergantung pada representasi Anda. Jika eksponen Anda sudah minimal sebelum mengurangi 1, metode ini akan memerlukan penanganan kasus khusus. Eksponen minimal pada xakan membuat x/2lebih kecil dari yang dapat diwakili (dengan asumsi mantissa diwakili dengan terkemuka 1).
    • Alih-alih mengurangi 1 dari eksponen x, xmantissa bergeser ke kanan dengan satu (dan tambahkan memimpin 1 implisit, jika ada).
    • Kurangi 1 dari eksponen y, jika tidak minimal. Jika minimal (y lebih besar dari x, karena mantissa), geser mantissa ke kanan dengan satu (tambahkan lead 1 implisit, jika ada).
    • Geser mantissa baru xke kanan sesuai dengan eksponen y.
    • Lakukan penambahan bilangan bulat pada mantissae, kecuali jika mantissa xtelah bergeser sepenuhnya. Jika kedua eksponen minimal, yang memimpin akan meluap, yang ok, karena yang meluap seharusnya menjadi yang terdepan lagi.
  • dan satu penambahan floating point.
    • Tidak dapat memikirkan kasus khusus apa pun di sini; kecuali untuk pembulatan, yang juga berlaku untuk pemindahan yang dijelaskan di atas.
Tidak ada Jawaban
sumber