26 views (last 30 days)

I am encountering a precision error with Matlab2020b, which I did not have in version 2016b.

I have 78-dimension vector x (attached). if I do the following, even though the result should be 0, I get a complex number as a result from acos calculation:

> y = x;

> acos(dot(x,y)/sqrt(sum(x.^2)*sum(y.^2)))

ans = 0.0000e+00 + 2.1073e-08i

In Matlab2016b, I know that using "norm" function caused a precision error and acos(dot(x,y)/(norm(x)*norm(y)) gave a complex number.

Back then, the use of sqrt(sum(x.^2)*sum(y.^2)) was a recommended method to avoid this issue. (as summarized in this page: https://stackoverflow.com/questions/36093673/why-do-i-get-a-complex-number-using-acos)

This method has been working fine in 2016b, but now with exactly the same code I have the complex number issue coming back in 2020b.

Was there a change in the precision of calculation in the newer version of matlab? If so, is there any good work around to avoid this issue?

Thanks,

Hiroyuki

Bruno Luong
on 15 Oct 2020

Edited: Bruno Luong
on 15 Oct 2020

This is a robust code.

theta = acos(max(min(dot(x,y)/sqrt(sum(x.^2)*sum(y.^2)),1),-1))

Note it returns 0 for x or y is 0. One might prefer NaN because correlation is undefined.

Jan
on 20 Oct 2020

The ACOS function is numerically instable at 0 and pi.

SUM is instable at all. A trivial example: sum([1, 1e17, -1]) .There are different approaches to increase the accuracy of the summation, see https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum

There is a similar approach for a stabilized DOT product, but the problem of ACOS will still exist. To determine the angle between two vectors, use a stable ATAN2 method, see https://www.mathworks.com/matlabcentral/answers/471918-angle-between-2-3d-straight-lines#answer_383392

Jan
on 21 Oct 2020

@Bruno: Sorry, my fault. I meant the "2*atan" approach, not "atan2":

x = rand(1, 78);

y = rand(1, 78);

angle1 = acos(dot(x, y) / (norm(x) * norm(y)))

angle2 = 2 * atan(norm(x * norm(y) - norm(x) * y) / ... % Typo fixed

norm(x * norm(y) + norm(x) * y))

Paul
on 22 Oct 2020

I think you have an error in angle2. Should it not be:

angle2 = 2*atan(norm(norm(x)*y - norm(y)*x)/norm(norm(x)*y + norm(y)*x))

Uday Pradhan
on 15 Oct 2020

Edited: Uday Pradhan
on 16 Oct 2020

Hi Hiroyuki,

If you check (in R2020b):

>> X = dot(x,y) - sqrt(sum(x.^2)*sum(y.^2))

ans =

1.776356839400250e-15

where as in R2016b, we get:

>> dot(x,y) - sqrt(sum(x.^2)*sum(y.^2))

ans =

0

Hence, in R2020b, we get:

>> acos(X)

ans =

0.000000000000000e+00 + 2.107342425544702e-08i

This is because the numerator dot(x,y) is "greater" than the denominator sqrt(sum(x.^2)*sum(y.^2)) albeit by a very small margin and hence the fraction X becomes greater than 1 and thus acos(X) gives complex value.

To avoid this my suggestion would be to establish a threshold precision to measure equality of two variables, for example you could have a check function so that if abs(x-y) < 1e-12 then x = y

function [a,b] = check(x,y)

if abs(x-y) < 1e-12

a = x;

b = a;

end

end

Now, you can do [a,b] = check(x,y) and then call acos(a/b). This will also help in any other function where numerical precision can cause problems.

Hope this helps!

Bruno Luong
on 18 Oct 2020

Something weird that I don't quite understand. In v2020b, the dot code is

...

if isreal(a) && isreal(b)

c = a'*b;

return

end

elseif ~isequal(size(a),size(b))

error(message('MATLAB:dot:InputSizeMismatch'));

end

c = sum(conj(a).*b);

...

Why calculate using a'*b for a, b reals and sum(conj(a).*b) for complex?

If I was TMW I would code a'*b for both cases. It might have something to do with interleave complex internal storage, but I fail to see the advantage of using sum(conj(a).*b)

Paul
on 18 Oct 2020

Paul
on 19 Oct 2020

If you look further down in dot.m (2019a) to the section used when dim is specified, you will see that there is a path to

c = sum(conj(a).*b,dim)

even if a and b are both vectors and are isreal. For example

dot(1:3,1:3,1)

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 1 Comment

## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/612051-calculation-precision-changed-in-2020b#comment_1061108

⋮## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/612051-calculation-precision-changed-in-2020b#comment_1061108

Sign in to comment.