Cara komputasi yang stabil secara numerik antara vektor

Saat menerapkan rumus klasik untuk sudut antara dua vektor:

α = \arccos \frac{v_{1} \cdot v_{2}}{‖ v_{1} ‖ ‖ v_{2} ‖}

$\alpha = \arccos \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{\|\mathbf{v_1}\| \|\mathbf{v_2}\|}$

orang menemukan bahwa, untuk sudut yang sangat kecil / akut, ada kehilangan presisi dan hasilnya tidak akurat. Seperti yang dijelaskan dalam jawaban Stack Overflow ini , salah satu solusinya adalah menggunakan arctangent sebagai gantinya:

α = \arctan 2 (‖ v_{1} \times v_{2} ‖, v_{1} \cdot v_{2})

$\alpha = \arctan2 \left(\|\mathbf{v_1} \times \mathbf{v_2}\|, \mathbf{v_1} \cdot \mathbf{v_2} \right)$

Dan ini memang memberikan hasil yang lebih baik. Namun, saya bertanya-tanya apakah ini akan memberikan hasil yang buruk untuk sudut yang sangat dekat dengan $\pi / 2$ . Apakah ini masalahnya? Jika demikian, apakah ada rumus untuk menghitung sudut secara akurat tanpa memeriksa toleransi di dalam ifcabang?

algorithms precision numerical-limitations astrojuanlu
sumber

Ini akan tergantung pada implementasi dari fungsi tangen invers dua parameter. Versi lambat dan stabil beralih secara kondisional antara bekerja dengan x / y dan y / x untuk mempertahankan presisi, sementara yang cepat hanya memasukkan hal-hal di kuadran kanan dan dengan demikian tidak ada yang lebih tepat daripada versi satu parameter.

origimbo

Anda harus mendefinisikan "kehilangan presisi": misalkan jawaban yang tepat adalah

α

$\alpha$ dan Anda mendapatkan

α + Δ

$\alpha + \Delta$ . Apakah Anda memerlukan

Δ ≪ α

$\Delta \ll \alpha$ atau

Δ ≪ π

$\Delta \ll \pi$ sudah cukup?

Stefano M

Dalam hal ini, jawaban yang benar adalah dan saya mendapat , keduanya .

α

$\alpha$

α \cdot 10^{8}

$\alpha \cdot 10^8$

≪ 1

$\ll 1$

astrojuanlu

Jawaban:

( Saya telah menguji pendekatan ini sebelumnya, dan saya ingat itu berhasil dengan benar, tetapi saya belum mengujinya secara khusus untuk pertanyaan ini. )

Sejauh yang saya tahu, keduanya dan dapat menderita pembatalan besar jika hampir paralel / tegak lurus — atan2 tidak dapat memberikan akurasi yang baik jika salah satu input dimatikan. $\|\mathbf{v}_1\times \mathbf{v}_2\|$ $\mathbf{v}_1\cdot \mathbf{v}_2$

Mulailah dengan merumuskan kembali masalah sebagai menemukan sudut segitiga dengan panjang sisi,dan(Ini semua dihitung secara akurat dalam aritmatika floating point). Ada varian rumus Heron yang terkenal karena Kahan ( Area salah hitung dan Sudut Segitiga yang Seperti Jarum ), yang memungkinkan Anda menghitung luas dan sudut (antara dan ) dari sebuah segitiga yang ditentukan oleh panjang sisinya, dan melakukannya secara stabil secara numerik. Karena pengurangan subproblem ini juga akurat, pendekatan ini harus bekerja untuk input yang sewenang-wenang. $a=|\mathbf{v}_1|$ $b=|\mathbf{v}_2|$ $c=|\mathbf{v}_1-\mathbf{v}_2|$ $a$ $b$

Mengutip dari makalah itu (lihat hal.3), dengan asumsi , Semua kurung di sini ditempatkan dengan hati-hati, dan mereka peduli; jika Anda menemukan diri Anda mengambil akar kuadrat dari angka negatif, panjang sisi input bukan panjang sisi segitiga. $a\geq b$

μ = {\begin{cases} c - (a - b), & if b \geq c \geq 0, \\ b - (a - c), & if c > b \geq 0, \\ invalid triangle, & otherwise \end{cases}

$\mu = \begin{cases} c-(a-b),&\text{if }b\geq c\geq 0,\\ b-(a-c),&\text{if }c>b\geq 0,\\ \text{invalid triangle},&\text{otherwise} \end{cases}$

a n g l e = 2 \arctan (\sqrt{\frac{((a - b) + c) μ}{(a + (b + c)) ((a - c) + b)}})

$\mathrm{angle} = 2\arctan\left( \sqrt{\frac{((a-b)+c)\mu}{(a+(b+c))((a-c)+b)}} \right)$

Ada penjelasan tentang bagaimana ini bekerja, termasuk contoh nilai yang gagal rumus lainnya, dalam makalah Kahan. Formula pertama Anda untuk adalah di halaman 4. $\alpha$ $C''$

Alasan utama saya menyarankan rumus Kahan's Heron adalah karena itu membuat primitif yang sangat bagus — banyak pertanyaan geometri planar yang berpotensi rumit dapat direduksi menjadi menemukan area / sudut segitiga sembarang, jadi jika Anda dapat mengurangi masalah Anda menjadi itu, ada formula stabil yang bagus untuk itu, dan tidak perlu untuk membuat sesuatu sendiri.

Sunting Mengikuti komentar Stefano, saya membuat plot kesalahan relatif untuk , ( kode ). Dua baris tersebut adalah kesalahan relatif untuk dan , sepanjang sumbu horizontal. Tampaknya itu berfungsi. $v_1=(1,0)$ $v_2=(\cos\theta, \sin\theta)$ $\theta=\epsilon$ $\theta=\pi/2-\epsilon$ $\epsilon$

Kirill
sumber

Terima kasih atas tautan dan jawabannya! Sayangnya formula kedua yang saya tulis tidak muncul di artikel. Di sisi lain, metode ini bisa sedikit rumit, karena memerlukan proyeksi dalam 2D.

astrojuanlu

@astrojuanlu Tidak ada proyeksi untuk 2d di sini: apa pun dua vektor 3d itu, mereka mendefinisikan satu segitiga (planar) di antara mereka — Anda hanya perlu tahu panjang sisinya.

Kirill

Anda benar, komentar saya tidak masuk akal. Saya berpikir dalam koordinat bukannya panjang. Terima kasih lagi!

astrojuanlu

@astrojuanlu Satu hal lagi yang ingin saya catat: tampaknya ada bukti formal bahwa rumus area akurat dalam Cara Menghitung Luas Segitiga: Kunjungan Resmi , Sylvie Boldo , menggunakan Flocq.

Kirill

Jawaban yang sangat bagus, tetapi saya berpendapat bahwa Anda selalu dapat menghitung secara akurat dalam aritmatika floating point. Bahkan jika maka pembatalan katastropik terjadi dalam menghitung komponen .

c

$c$

c < ϵ \cdot min (a, b)

$c < \epsilon \cdot \min(a,b)$

(v_{1} - v_{2})

$(\mathbf{v}_1-\mathbf{v}_2)$

Stefano M

Jawaban efisien untuk pertanyaan ini, tidak terlalu mengejutkan, dalam catatan lain oleh Velvel Kahan :

α = 2 \arctan (‖ \frac{v_{1}}{‖ v_{1} ‖} + \frac{v_{2}}{‖ v_{2} ‖} ‖, ‖ \frac{v_{1}}{‖ v_{1} ‖} - \frac{v_{2}}{‖ v_{2} ‖} ‖)

$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right)$

di mana saya menggunakan sebagai sudut yang dibuat oleh dengan sumbu horizontal. (Anda mungkin harus membalik urutan argumen dalam beberapa bahasa.) $\arctan(x,y)$ $(x,y)$

(Saya memberikan demonstrasi Mathematica dari formula Kahan di sini .)

JM
sumber

Maksudmu ?

\arctan 2

$\arctan2$

astrojuanlu

Saya terbiasa menggambarkan arctangent dua argumen sebagai , ya. Dalam bahasa seperti FORTRAN, padanannya adalah .

\arctan (x, y)

$\arctan(x,y)$ ATAN2(Y, X)