[Numpy-discussion] Arcos returns nan

Charles R Harris charlesr.harris at gmail.com
Fri Mar 28 12:05:56 EDT 2008


On Fri, Mar 28, 2008 at 3:20 AM, João Quinta da Fonseca <
joao.q.fonseca at gmail.com> wrote:

> I have a function that returns the dot product of two unit vectors.
> When I try to use arcos on the values returned I sometimes get the
> warning:
> "Warning: invalid value encountered in arccos", and the angle
> returned is nan. I found out that this happens for essentially co-
> linear vectors, for which the dot product function returns 1.0.  This
> looks like 1.0 no matter how I print it but if I do: >>N.dot(a,b)>1,
> I get: >>True.


Then the dot product is > 1. If you print out the number with around 20
significant digits you will see that. The differences between Matlab,
Fortran, and numpy may be due to compiler flags, the compile, and your
hardware. The internal registers in the floating point unit can have extra
precision, so how those registers are used can effect the low order bits of
the result. We really need to have your input vectors to see what is going
on here.


> Now I guess this arises because of the way computers store floats but
> shouldn't numpy take care of this somehow? Is it a bug? I don't seem
> to have this problem with Matlab or Fortran.
>

You need to check for the correct domain, as using the arccos for this when
the vectors are nearly colinear is going to be tricky and sensitive to
roundoff just from the nature of the problem. Arccos is also not very
accurate for values near in any case +/- 1 because the extremum of the cos
function has those values. If you are working in 2 or three dimensions you
could try the cross product, or more generally, you can normalize the
vectors and look at the difference vector, which will solve the nan problem,
but it won't give you better precision.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080328/0b6ef01d/attachment.html>


More information about the NumPy-Discussion mailing list