DDtheta: numerical stability of small angular separation in float32

[JonLoveday]

General information

• Corrfunc version: 2.4.0
• platform: MacOS
• installation method (pip/source/other?): pip

Issue description

When I perform pair counts on a random distribution with float32 precision using Corrfunc.mocks.DDtheta_mocks with 20 log bins between 0.01 and 10 degrees, I always get zero pairs in the second bin. Using float64 alleviates this problem. I'm not sure if this is a bug as such, but users should be strongly advised to only use float64 coordinates.

Expected behavior

Would not expect significantly different results betweein flloat32 and float64 coords for angular precision > 0.01 deg.

Actual behavior

Zero pair counts in second bin using float32.

What have you tried so far?

Using float64 precision avoids this problem. In fact, pair counts in all bins change.

Minimal failing example

import Corrfunc
import numpy as np
from numpy.random import default_rng
rng = default_rng()
def corrfunc_test(nran=10000, tmin=0.01, tmax=10, nbins=20):
    """Random-random pair count test."""

    bins = np.logspace(np.log10(tmin), np.log10(tmax), nbins + 1)
    ra = 20*rng.random(nran)
    dec = 20*rng.random(nran)
    counts = Corrfunc.mocks.DDtheta_mocks(1, 1, bins, ra, dec)
    print(counts)

    counts = Corrfunc.mocks.DDtheta_mocks(1, 1, bins, ra.astype('float32'), dec.astype('float32'))
    print(counts)

[lgarrison]

At first I thought this might be due to acos precision, but now I don't think it is.

@manodeep I think the method we're using to compute delta-theta is not numerically stable. We're doing

DOUBLE costheta = (one - half*sqr_chord_sep);

where sqr_chord_sep will be approximately (0.02/180*np.pi)**2 = 1.2e-7 for an angular separation of 0.02 deg. Then the costheta subtraction underflows to 1, unsurprisingly. It seems to me we should be able to do better than this.

Perhaps we could switch to theta = 2*arcsin(C/2) at some small angle; the expression is exact, but slower because it involves both arcsin and sqrt.

Or perhaps we could switch to the small angle approximation theta = C. For 0.02 deg, the error is only 2e-12. The next term in the expansion is C^3/24, which is also cheap to compute once we have C. It should be possible to compute the angle for which the error due to 1-C^2 is greater than the error due to the small angle approximation, and use that as the switchover point.

Or perhaps there's another trig trick that avoids 1-C^2?

[manodeep]

Thanks @JonLoveday for reporting the issue! It's not entirely unexpected - computing with floating points numbers can cause such seemingly odd behaviour. Corrfunc was not designed to be resilient against such numerical issues - for example, rpavg/weightavg values can differ between float and double precision.

What I would also be curious to know is if a brute-force python computation with float-precision produced results more consistent with double-precision.

@lgarrison Yeah $DD(\theta)$ is definitely the most susceptible to such floating point issues. Have not really thought about or investigated any numerical tricks to minimise such "catastrophic" errors. I have a vague memory that $DD(\theta)$ can compute theta ranges from [0, 180] - any numerical tricks we use should not reduce the max-theta range allowed. May be we could just use a python script to write out the absolute and relative errors for all "small" angles and then decide the switchover point

[JonLoveday]

Thanks @lgarrison and @manodeep! In the past, I have always used theta = 2*arcsin(C/2). It just occurred to me that you can avoid both the arcsin and sqrt calculations by converting the bins from theta -> 4 sin^2(theta/2).

[manodeep]

Currently Corrfunc can compute angular separations between [0, 180] - i.e., $-1.0 <= cos(\theta) <= 1.0$. If we take the squared route, then we will be limited to [0, 90] - which may not suit what others want. However, we would likely only use the squared approx. near zero separation anyway - so perhaps the entire angular range would not be affected.

[JonLoveday]

sin^2(theta/2) is a unique mapping of theta over [0, 180]

[manodeep]

Right! $\theta/2$ and not $\theta$ 🤦🏾

[lgarrison]

Thanks for the great suggestion, @JonLoveday! I made a proof-of-concept in #297. The biggest hurdle remaining is probably deciding how to compute theta for purposes of returning the average binned theta. We could bite the bullet and do 2*arcsin(sqrt(C^2)/2)... but for most angles arccos(1 - C^2/2) is probably accurate enough, and theta = C might suffice for small angles (maybe with an additional C^3/24 term).

On the other hand, whatever we choose has to be implemented in the SIMD kernels, which makes a theta computation with two branches less appealing.

I've updated the title to reflect the underlying issue.

[JonLoveday]

Thanks for the great suggestion, @JonLoveday! I made a proof-of-concept in #297. The biggest hurdle remaining is probably deciding how to compute theta for purposes of returning the average binned theta. We could bite the bullet and do 2*arcsin(sqrt(C^2)/2)... but for most angles arccos(1 - C^2/2) is probably accurate enough, and theta = C might suffice for small angles (maybe with an additional C^3/24 term).

I think for most practical purposes, it would be good enough to calculate the average binned C^2 value, and then transform to an average theta.

[manodeep]

While looking at #301, I realised that we could easily recast the current calculation to be:
costheta := x1*x2 + y1*y2 + z1*z2 and remove the entire chord_sep calculation. That should be better behaved and much easier to implement within the existing framework

[lgarrison]

I think this issue is less about avoiding the chord_sep calculation and more about choice of variables. If cos(theta) is the variable, no matter how you compute it, it's going to have very poor precision for small angles. A dot product wouldn't fix the original issue, for example—the cosine of an angle less than ~0.02 degrees is indistinguishable from 1 in float32, regardless of how you get there. So I don't think we'd see much gain from a dot-product approach, but I'd be happy to be proven wrong!

[manodeep]

I thought that the issue was arising because 1 - 0.5*C^2 for small C was numerically unstable, and switching to the dot product would add terms that were all of similar magnitude, and hence better-behaved numerically.

@lgarrison I can confirm that your assertion is correct - even with the dot-product method (in #303), the second bin still contains 0 elements for a float32 calculation. Now I understand what is going on - the difference between the two bin-edges for the second bin is 3e-8 and is smaller than float-eps (few*10^-7) - therefore, the bin can not be populated with a float32. Maybe that should be the criteria to switch to a different variable or even check that condition is satisfied or even "refuse" to compute in float32?

[manodeep]

Modifying the reproducer slightly to check the bins:

import numpy as np
def corrfunc_test(nran=10000, tmin=0.01, tmax=10, nbins=20):
    """Random-random pair count test."""

    bins = np.geomspace(tmin, tmax, nbins + 1)
    for typ in [np.float32, np.float64]:
        cosbins = np.cos(bins*np.pi/180., dtype=typ)
        diff = np.diff(cosbins)
        print("float type = {} max. diff = {}".format(typ, max(diff)))

produces

float type = <class 'numpy.float32'> max. diff = 0.0
float type = <class 'numpy.float64'> max. diff = -1.5158711841323225e-08