Calculation precision changed in 2020b?

3 visualizzazioni (ultimi 30 giorni)
Hiroyuki Kato
Hiroyuki Kato il 13 Ott 2020
Commentato: Jan il 23 Ott 2020
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
  1 Commento
Bruno Luong
Bruno Luong il 15 Ott 2020
Modificato: Bruno Luong il 15 Ott 2020
"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."
Pure luck. None of the observation has rigorous justification.

Accedi per commentare.

Risposta accettata

Bruno Luong
Bruno Luong il 15 Ott 2020
Modificato: Bruno Luong il 15 Ott 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.

Più risposte (2)

Jan
Jan il 20 Ott 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
  4 Commenti
Paul
Paul il 22 Ott 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))
Jan
Jan il 23 Ott 2020
@Paul: Yes, there was a typo.

Accedi per commentare.


Uday Pradhan
Uday Pradhan il 15 Ott 2020
Modificato: Uday Pradhan il 16 Ott 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.
Another workaround can be found in this link : Determine the angle between two vectors.
Hope this helps!
  10 Commenti
Paul
Paul il 18 Ott 2020
Modificato: Paul il 18 Ott 2020
I will assume that the decision to change the implementation and yield a different answer was not undertaken lightly. I checked the release notes and did not see this change to dot, though the sum function did change in 2020B.
Paul
Paul il 19 Ott 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)

Accedi per commentare.

Categorie

Scopri di più su MATLAB in Help Center e File Exchange

Tag

Prodotti


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by