3
$\begingroup$

I am calculating the angle between vectors $$cos(\theta{_a}{_b}) = \frac{a^T b}{||a|| ||b||}$$

In mathematica I am coding

ArcCos[N[Transpose[{{10},{12},{4}}].{{2.5},{3},{1}}/(Norm[{{10},{12},{4}}] * Norm[{{2.5},{3},{1}}])][[1]][[1]]]

The dot product on top and and product of the norms on the bottom both equal 65. which means I should be calculating the arccos of one. The angle should be zero.

The answer I get is $$0 + 2.10734 x 10 ^{-8} i $$

If I use a fraction instead of decimal ...

ArcCos[Transpose[{{10},{12},{4}}].{{5/2},{3},{1}}/(Norm[{{10},{12},{4}}]*Norm[{{5/2},{3},{1}}])]

I get a zero.

What do I need to do to get a clean zero with decimals in my vectors?

Improve this question
$\endgroup$
2
  • $\begingroup$ Welcome to the Mathematica Stack Exchange. You can use the two argument form of Chop: ArcCos[N[Transpose[{{10}, {12}, {4}}] . {{2.5}, {3}, \ {1}}/(Norm[{{10}, {12}, {4}}]*Norm[{{2.5}, {3}, {1}}])][[1]][[1]]] // Chop[#, 10^-7] & $\endgroup$ Commented Feb 25, 2023 at 14:38
  • 1
    $\begingroup$ Have you seen VectorAngle? Alternatively, use a numerically stable formula such as With[{a = Flatten@N@{{10}, {12}, {4}}, b = Flatten@N@{{2.5}, {3}, {1}}}, With[{n1 = Norm[a], n2 = Norm[b]}, (*vector angle formula by Velvel Kahan*) 2 ArcTan[Norm[a n2 + n1 b], Norm[a n2 - n1 b]] ]] from here. $\endgroup$ Commented Feb 25, 2023 at 17:17

4 Answers 4

Reset to default
4
$\begingroup$

First note, the expression:

Transpose[{{10}, {12}, {4}}] . {{2.5}, {3}, {1}}

does not make sense. You are trying to mimic column and row vectors, but MMA does not work with row and column vectors. Instead the dot product sums over the rightmost index of the left object and the leftmost index of the right object. All in all, the above expression in MMA is written simply by:

{10, 12, 4} . {2.5, 3, 1}

Your expression in correct MMA notation is then:

ArcCos[{10, 12, 4} . {2.5, 3, 1}/(Norm[{10, 12, 4}] Norm[{2.5, 3, 1}])]

what evaluates to zero.

On the other hand, Look what you are feeding to ArcCos:

N[Transpose[{{10}, {12}, {4}}] . {{2.5}, {3}, {1}}/(Norm[{{10}, {12}, \
{4}}]*Norm[{{2.5}, {3}, {1}}])][[1]][[1]] // FullForm

1.0000000000000002`

An argument larger than 1 will result in a complex ArcCos. Also note that the "N" is not needed as you already have machine numbers.

To fix this you can use "Chop" to get rid of spurious. But it is better to use correct syntax.

Improve this answer
$\endgroup$
1
  • $\begingroup$ One should note that a . b/(Norm[a] Norm[b]) is not guaranteed to be ≤ 1, if computed with machine precision. Example: a = {4.1, 5.34, 4.1}; b = (1 + $MachineEpsilon) a;. It does work on the OP's example, tho. $\endgroup$ Commented Feb 25, 2023 at 17:39
2
$\begingroup$

The numerical problems are two-fold and are independent of syntax and the OP's use of column matrices as vectors. The most serious problem is that round-off error leads to a range of $$ -1-\delta \le {a^T b\over\|a\|\,\|b\|} \le 1+\delta \,,$$ where $\delta$ is of the order of a small multiple of machine epsilon. And for values outside of $\pm1$, the arc cosine is imaginary. The other problem is that $\arccos$ is singular at $\pm1$,and an error of $\delta$ is propagated to an error of $\approx \sqrt{2\delta}$. At double precision, a round-off error in the 16th digit becomes an error in the 8th.

Here's a particularly bad example of a vector $a$ such that the "angle" between it and a positive scalar multiple $b=k\,a$ of itself is imaginary more than 75% of the time:

With[{a = {1.075`, 1.175`, 1.025`}},
  Table[With[{b = (1 + k*$MachineEpsilon) a},
      ArcCos[a . b/Norm[a]/Norm[b]]],
     {k, 1000}] // Im // Unitize
  ] // Total

(*  779  *)

A much better way is to use a numerically stable algorithm. Currently (V13.2), VectorAngle[] seems to do so. Tests show it does not use exactly the formula below, which was once the suggested replacement for VectorAngle[]. See Complex result for Real vectors in VectorAngle or Fast Spherical Linear Interpolation of list of quaternions Also, it may be adapted to a compilable formula by replacing Norm[v] by Sqrt[v.v].

The formula below computes the angle $\theta$ from the diagonals of a rhombus spanned by the normalizations of the vectors $a$ and $b$:

The use of the formula in numerical applications is credited to Velvel Kahan:

2 ArcTan[
 Norm[Norm[a] b + Norm[b] a],
 Norm[Norm[a] b - Norm[b] a]]

You cannot get a complex results since Norm[] is nonnegative. Further the error propagated is proportional to the round-off error. The value of the formula differs from that of VectorAngle[] occasionally, up to an error in the last digit of precision.

Improve this answer
$\endgroup$
1
$\begingroup$

This is not an answer, but too big for a comment. What you found is very annoying, and I can provide no solution. But maybe the following code is of some help.

a = {10, 12, 4};
b = {2.5, 3, 1};

The dot product

a.b

returns 2*Sqrt[65.]. The ratio

a.b/(Norm[a]*Norm[b])

returns 1., which is 1.0 to machine precision.

Checking this with Equal

a.b/(Norm[a]*Norm[b]) == 1.0

returns True, and checking this with SameQ

a.b/(Norm[a]*Norm[b]) === 1.0

also returns True. But

ArcCos[1.0]

returns 0., while

ArcCos[a.b/(Norm[a]*Norm[b])]

returns your answer with a small imaginary component. This imaginary component is even too big to get chopped off with the default call of Chop. Very annoying.

Improve this answer
$\endgroup$
1
$\begingroup$

This problem does not occur when normalizing the vectors with Normalize:

a = {10, 12, 4};
b = {2.5, 3, 1};

ArcCos[Normalize[a] . Normalize[b]]
(*    0.    *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Alas: SeedRandom[1]; aa = RandomReal[{-1, 1}, 3]; bb = (1 + 369*$MachineEpsilon) aa; ArcCos[Normalize[aa] . Normalize[bb]] returns Sqrt[2 $MachineEpsilon] I. But it's much better than the non-normalized formula. $\endgroup$ Commented Feb 26, 2023 at 20:46
  • $\begingroup$ @MichaelE2 Ah yes, very true! A minimal example would be a = {1, 1, 1} - $MachineEpsilon, for which Normalize[a] . Normalize[a] gives 1 + $MachineEpsilon. $\endgroup$ Commented Feb 27, 2023 at 9:53
  • $\begingroup$ It's really about getting the IEEE 754 rounding effects correct, which seems like an impossible task to get right in all cases. In the above example, $1/\sqrt{3}$ is represented as the double-precision number $\frac{5200308914369309}{9007199254740992}$ (you can see this with SetPrecision[N[1/Sqrt[3]], ∞]), which is less accurate than $\frac{5200308914369308}{9007199254740992}$ (decreasing the numerator by 1). A faulty rounding operation somewhere? $2^{53}/\sqrt{3}=5\,200\,308\,914\,369\,308.3...$ after all. $\endgroup$ Commented Feb 27, 2023 at 11:26
  • $\begingroup$ Maybe it's the difference between Sqrt[1/3.] and 1/Sqrt[3.]. $\endgroup$ Commented Feb 27, 2023 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.