Forma numéricamente estable de calcular ángulos entre vectores

14

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

α=arccosv1v2v1v2

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:

α=arctan2(v1×v2,v1v2)

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?π/2if

astrojuanlu
fuente
1
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? αα+ΔΔαΔπ
Stefano M
En este caso, la respuesta correcta fue y obtuve , ambos . αα1081
astrojuanlu

Respuestas:

18

( 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.v1×v2v1v2

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=|v1|b=|v2|c=|v1v2|ab

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.ab

μ={c(ab),if bc0,b(ac),if c>b0,invalid triangle,otherwise
angle=2arctan(((ab)+c)μ(a+(b+c))((ac)+b))

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.α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. v1=(1,0)v2=(cosθ,sinθ)θ=ϵθ=π/2ϵϵingrese la descripción de la imagen aquí

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
2
@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
2
@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 . cc<ϵmin(a,b)(v1v2)
Stefano M
7

La respuesta eficiente a esta pregunta es, no sorprendentemente, en otra nota de Velvel Kahan :

α=2arctan(v1v1+v2v2,v1v1v2v2)

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 ? arctan2
astrojuanlu
1
Estoy acostumbrado a representar el arcotangente de dos argumentos como , sí. En un lenguaje como FORTRAN, el equivalente sería . arctan(x,y)ATAN2(Y, X)
JM