Forma numéricamente estable de calcular ángulos entre vectores

Al aplicar la fórmula clásica para el ángulo entre dos vectores:

α = \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}\|}$

se encuentra que, para ángulos muy pequeños / agudos, hay una pérdida de precisión y el resultado no es exacto. Como se explica en esta respuesta de Stack Overflow , una solución es usar el arcotangente en su lugar:

α = \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)$

Y esto de hecho da mejores resultados. Sin embargo, me pregunto si esto daría malos resultados para ángulos muy cercanos a . Es el caso? Si es así, ¿hay alguna fórmula para calcular ángulos con precisión sin verificar la tolerancia dentro de una rama? $\pi / 2$ if

algorithms precision numerical-limitations astrojuanlu
fuente

Esto dependerá de la implementación de la función de tangente inversa de dos parámetros. Las versiones lentas y estables cambian condicionalmente entre trabajar con x / y e y / x para mantener la precisión, mientras que las rápidas simplemente pegan las cosas en el cuadrante derecho y, por lo tanto, no son más precisas que la versión de un parámetro.

origimbo

Debe definir "pérdida de precisión": suponga que la respuesta correcta es y obtiene en su lugar . ¿Necesita o es suficiente?

α

$\alpha$

α + Δ

$\alpha + \Delta$

Δ ≪ α

$\Delta \ll \alpha$

Δ ≪ π

$\Delta \ll \pi$

Stefano M

En este caso, la respuesta correcta fue y obtuve , ambos .

α

$\alpha$

α \cdot 10^{8}

$\alpha \cdot 10^8$

≪ 1

$\ll 1$

astrojuanlu

Respuestas:

( He probado este enfoque antes, y recuerdo que funcionó correctamente, pero no lo he probado específicamente para esta pregunta ) .

Por lo que puedo decir, ambos y puede sufrir una cancelación catastrófica si son casi paralelas / perpendiculares; atan2 no puede darle una buena precisión si alguna de las entradas está desactivada. $\|\mathbf{v}_1\times \mathbf{v}_2\|$ $\mathbf{v}_1\cdot \mathbf{v}_2$

Comience reformulando el problema como encontrar el ángulo de un triángulo con longitudes laterales,y(todos estos se calculan con precisión en aritmética de coma flotante). Hay una variante bien conocida de la fórmula de Heron debido a Kahan ( Area calcular mal y ángulos de una aguja similar a Triangle ), que le permite calcular el área y el ángulo (entre y ) de un triángulo especificado por sus longitudes de los lados, y hacerlo numéricamente estable. Debido a que la reducción a este subproblema también es precisa, este enfoque debería funcionar para entradas arbitrarias. $a=|\mathbf{v}_1|$ $b=|\mathbf{v}_2|$ $c=|\mathbf{v}_1-\mathbf{v}_2|$ $a$ $b$

Citando de ese documento (ver p.3), suponiendo , Todos los paréntesis aquí se colocan cuidadosamente, y son importantes; si te encuentras tomando la raíz cuadrada de un número negativo, las longitudes laterales de entrada no son las longitudes laterales de un triángulo. $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)$

Hay una explicación de cómo funciona esto, incluidos ejemplos de valores para los que fallan otras fórmulas, en el artículo de Kahan. Su primera fórmula para es en la página 4. $\alpha$ $C''$

La razón principal por la que sugiero la fórmula de Kahan Heron es porque es una primitiva muy agradable: muchas preguntas de geometría plana potencialmente complicadas se pueden reducir a encontrar el área / ángulo de un triángulo arbitrario, por lo que si puede reducir su problema a eso, hay es una fórmula estable y agradable, y no hay necesidad de inventar algo por su cuenta.

Editar Siguiendo el comentario de Stefano, hice un diagrama de error relativo para , ( código ). Las dos líneas son los errores relativos para y , van a lo largo del eje horizontal. Parece que funciona. $v_1=(1,0)$ $v_2=(\cos\theta, \sin\theta)$ $\theta=\epsilon$ $\theta=\pi/2-\epsilon$ $\epsilon$

Kirill
fuente

Gracias por el enlace y la respuesta! Lamentablemente, la segunda fórmula que escribí no aparece en el artículo. Por otro lado, este método puede ser un poco complejo, ya que requiere proyección en 2D.

astrojuanlu

@astrojuanlu Aquí no hay proyección para 2d: sean cuales sean los dos vectores 3d, definen un único triángulo (plano) entre ellos; solo necesita saber sus longitudes laterales.

Kirill el

Tienes razón, mi comentario no tiene sentido. Estaba pensando en coordenadas en lugar de longitudes. ¡Gracias de nuevo!

astrojuanlu

@astrojuanlu Una cosa más que quiero señalar: parece que hay una prueba formal de que la fórmula del área es precisa en Cómo calcular el área de un triángulo: una revisión formal , Sylvie Boldo , usando Flocq.

Kirill

Excelente respuesta, pero discrepo de que siempre se puede calcular con precisión en aritmética de coma flotante. De hecho, si se producen cancelaciones catastróficas al calcular los componentes de .

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

La respuesta eficiente a esta pregunta es, no sorprendentemente, en otra nota de 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)$

donde uso como el ángulo formado por con el eje horizontal. (Es posible que deba cambiar el orden de los argumentos en algunos idiomas). $\arctan(x,y)$ $(x,y)$

(Di una demostración de Mathematica de la fórmula de Kahan aquí ).

JM
fuente

¿Te refieres a ?

\arctan 2

$\arctan2$

astrojuanlu

Estoy acostumbrado a representar el arcotangente de dos argumentos como , sí. En un lenguaje como FORTRAN, el equivalente sería .

\arctan (x, y)

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