Calculation precision changed in 2020b?

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]

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

Answer 1

[Bruno Luong]

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.

Answer 2

[Jan]

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

[Bruno Luong]

But the atan2 only applies for 2D and 3D (real) vectors. Here OP's example is in R^78, how do you apply in this case?

In defend of acos, it depends on what one do, sometime acos precision is sufficient with faster, simpler code. It's directly linked to scalar product and cauchy schwarz inequality.

It's good to be awared about precision limitation of each method, but not constantly taking atan2 method as a fixated objective.

[Jan]

@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))

See: http://people.eecs.berkeley.edu/~wkahan/Mindless.pdf and https://possiblywrong.wordpress.com/2020/07/17/computing-the-angle-between-two-vectors/

[Paul]

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]

@Paul: Yes, there was a typo.

Answer 3

[Uday Pradhan]

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!

[Paul]

Which function(s) used in those computations changed between those two Matlab versions?

[Uday Pradhan]

Hi Paul,

The method of handling real and complex arguments to "dot" has been updated in the current version of MATLAB which is resulting in this difference.

[Bruno Luong]

Sorry but your explanation is not accurate: in the stackoverflow page the dot result is exactly -3 regardless the version, both vectors has integer components.

The difference is surely from the cauculation of the denominator using NORM and SQRT(SUM...)

In anycase one cannot rely on math inequality that can only hold in ideal word, here we deal with finite precision floating point.

[Uday Pradhan]

>> dot(x,x)  %in MATLAB R2020b

ans =

  15.997157281555475
  
>> dot(x,x)     %in MATLAB R2016b

ans =

  15.997157281555474 

This is just an example, as I said there's a slight modification to handle real vectors. We get the same result for norm(x) and the sqrt (sum(...)) in both versions. I have also added a second workaround in my answer.

[Bruno Luong]

That just shows that the difference of results can be sourced from everywhere when a sum of three terms is involved. Program a code based on "it works" on a single example without justification or understading is insanely dangerous.

The robust approach is truncated the acos argument in [-1,1] as your "workaround" (to me it's not a workaround, it's the only right way to do it)

[Paul]

I don’t have either version in question here, so I’d still like to know: did the m-code in dot change, or did the implementation of the function(s) called by dot change, or both?

[Bruno Luong]

dot(x,y) simply calls

x'*y

in 2020

and

sum(conj(x).*y)

in 2016b

[Bruno Luong]

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]

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]

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)