Handling Numerical Instability in Estimating Angle Between 2D Points

Hi! I would like to estimate the angle between points using the following code.

x = randn(10000, 1); y = randn(10000, 1);
angle_mat = zeros(length(x));
for i = 1:length(x)
    for j = 1:length(x)
        if i ~= j
            angle_mat(i,j) = atan2(y(j)-y(i), x(j)-x(i));
        else
            angle_mat(i,j) = NaN;
        end
    end
end

I realized, however, that, due to numerical instability near horizontal/vertical slopes, the angle estimate near 0/+-pi became weird. Does anyone know how to acount for this numerical instability?

[Walter Roberson]

By the way, the code can be written without loops:

 x = randn(10000, 1); y = randn(10000, 1);
 angle_mat = atan2(y - y', x - x.');
 angle_mat(1:size(x,1)+1:end) = NaN;

[Walter Roberson]

I do not observe any particular problem ?

x = linspace(1-1e-5,1+1e-5,100);
y = x;
angle_mat = atan2(y - y', x - x.');
angle_mat(1:size(angle_mat,1)+1:end) = NaN;
 
 
 surf(x, y, angle_mat, 'edgecolor', 'none')

●

[Leo]

Hi, thank you! You are right; the simulated dataset I provided does not demonstrate the problem I posed. The problem, however, exists in my dataset, where there seems to be an asymptote-like behaviour near 0 and +- pi. Originally, I thought it was a numerical instability issue like explained in here, for example: https://ch.mathworks.com/matlabcentral/answers/78931-calculate-the-angle-between-multiple-points. But, perhaps it might just be a product of my dataset?

[William Rose]

@Leo, is the plot above a histogram of your data?

[edit: I now wonder if you are looking at little bumps on the histogram at 0 and +pi.]

Can you provide the code and the dataset that gave the plot above?

The big peaks in your plot are not at +-pi/2 and are not at pi. There is a small bump around 0, but would need to see an enlarged plot to see it this is really at 0. atan() and atan2() are well behaved at theta=0. There is an even smaller bump around theta=+pi. You'd have to look at the detials to understand if this bump is statistical noise or not. Could it be a rounding effect? For example, if x and y values are rounded to 3 significant figures, maybe the rounding, near some edges, does something funny to the distribution, causing slightly more points to round to the +pi side than to the -pi side. (The first bin on the -pi side is a bit low compared to its neighbors.) You could see if it is possible to obtain the data with more significant figures, and see if these bumps change.

Answer 1

[William Rose]

@Leo,

[moved my comment to be an answer, which is what I intended]

You are defining numerical instability as a discontinuity. The plot below shows the angle from the origin to a point that is orbiting the origin CCW at radius 1, starting at (x,y)=(1,0). atan2() has one discontinuity, which occurs when the orbiting point crosses the minus-x axis ("due west" of the origin). atan() has two discontinuities, at pi/2 and -pi/2. Use unwrap() to remove the discontinuity in atan2().

theta=0:pi/50:2*pi;
x=cos(theta); y=sin(theta);
plot(theta,atan(y./x),'-ro',theta,atan2(y,x),'-gx',theta,unwrap(atan2(y,x)),'-b+')
legend('atan','atan2','unwrap(atan2)')

●

[Leo]

Hi William, thank you for your answer! I am not sure if I am understanding your solution correctly, but I believe unwrap() here is not the answer I am seeking for. Due to the nature of the data, I cannot share it, but I can try to explain my algorithm the best I can.

amat = atan2(y - y', x - x.');
amat(1:sum(idx==1)+1:end) = NaN;

Essentially, by running, this code I am hoping to find the pairwise angle (a proxy for slope) between all my data points given by x and y. I hypothesise that there are two clusters with two different slopes, as you can see from the two peaks in my histogram. There, however, is an annoying peak at 0 and +/-pi, which I believe to due to the numerical instability when trying to calculate atan2() from a near infinite or near 0 slope. Here you can see it more emphasized in the same histogram but wrapped between [-pi/2 pi/2] (therefore the peaks are now at 0 and +/- pi/2).

I was wondering if there is a better code to prevent the error I am facing. Or do you know if I am mistakening my issue with another problem?

[William Rose]

@Leo,

Your angle histograms indicate more angles very close to 0 and +-pi/2 than would expect by chance. This suggests something is causing points to be horizontally aligned (same or almost same y-values) more often than expected. I would check the experimental procedures and data processing, including how the x and y values are stored and saved in files. Look for something such as rounding or quantization error that could be causing a slight enrichment in the number of identical y-values. I don;t think it is related to atan2(). In the example below, I make 6 plots of atan2() values, for 6 different sets of 200 points randomly scattered on the unit square. (I made six just to illustrate that the slight peaks and valleys across the histograms are different each time, i.e. random.)

figure;
for i=1:6
    x=rand(200,1); y=rand(200,1);
    angleMat=atan2(y-y',x-x');
    angleMat(1:size(x,1)+1:end) = NaN;
    subplot(2,3,i); histogram(angleMat)
end

There is no sign of excess values at 0 or at +pi/2 or at -pi/2. The outermost bins are a bit lower than the others, in each case. This is because the default bin edges extend a bit past +-pi/2, so the edge bins don't get as many hits as the others. This is expected.

Repeat, but this time, round the y-values to 1 digit after the decimal point.

figure;
for i=1:6
    x=rand(200,1); y=round(rand(200,1),1);
    angleMat=atan2(y-y',x-x');
    angleMat(1:size(x,1)+1:end) = NaN;
    subplot(2,3,i); histogram(angleMat)
end

●

It's not exactly like the spikes in your histogram, but is similar in some respects.

I would look at the histograms of x-x' and of y-y' for your data, to see if they have any unexpected features.