Fix consistentcy of histc on CPU and CUDA

[Aidyn-A]

Fixes #87657

The main reason why histc returns slightly different outputs is the difference on how bin position is calculate.
The CPU calculates it as:

pytorch/aten/src/ATen/native/cpu/HistogramKernel.cpp

Lines 168 to 170 in 449778a

	pos = static_cast<int64_t>((elt - leftmost_edge[dim])
	/ (rightmost_edge[dim] - leftmost_edge[dim])
	* (num_bin_edges[dim] - 1));

which is basically (i - a) / (b - a) * N, while cuda code

pytorch/aten/src/ATen/native/cuda/SummaryOps.cu

Line 41 in 449778a

IndexType bin = (int)(((bVal - minvalue)) * nbins / (maxvalue - minvalue));

which is (i - a) * N / (b - a).

For some cases like in #87657 the order of arithmetic operations matters due to the floating point round-off.

---

Not sure where would be the most appropriate place to put the unit test. Hope test_reductions::test_histc will do.

cc @VitalyFedyunin @jgong5 @mingfeima @XiaobingSuper @sanchitintel @ashokei @jingxu10

[pytorch-bot]

🔗 Helpful Links

🧪 See artifacts and rendered test results at hud.pytorch.org/pr/87832

• 📄 Preview Python docs built from this PR
• 📄 Preview C++ docs built from this PR
• ❓ Need help or want to give feedback on the CI? Visit our office hours

Note: Links to docs will display an error until the docs builds have been completed.

⏳ No Failures, 1 Pending

As of commit e187347:
💚 Looks good so far! There are no failures yet. 💚

This comment was automatically generated by Dr. CI and updates every 15 minutes.

[soumith]

great, if you fix lint we can merge it.

[ngimel]

would test_against_np catch it instead?

[Aidyn-A]

Unfortunately, there is a discrepancy with numpy for this input tensor.
The code

import torch
import numpy as np

def fn(x):
    return torch.histc(x, bins=10, min=0, max=0.99)

def fn_np(x):
    return torch.from_numpy(np.histogram(x.numpy(), bins=10, range=(0, 0.99))[0])


x = torch.linspace(0, 0.99, 1001, dtype=torch.float32)
torch.histc(x, bins=10, min=0, max=0.99)

print(fn(x))
print(fn(x.cuda()))
print(fn_np(x))

Produces:

tensor([101.,  99., 101.,  99., 101.,  99., 100., 100., 100., 101.])
tensor([101.,  99., 101.,  99., 101.,  99., 100., 100., 100., 101.], device='cuda:0')
tensor([101,  99, 100, 100, 100, 100, 100, 100, 100, 101])

[mruberry]

Is NumPy's histogram more accurate? Making histc consistent on CPU and CUDA is nice, but is the histc algorithm we're using actually correct?

[Aidyn-A]

NumPy might be more accurate because it uses direct comparison https://github.com/numpy/numpy/blob/13d55a3c2f016a58a6e9d6b8086f338e07c7478f/numpy/lib/histograms.py#L862-L866 to find out if the value is within the range or not. However, I noticed that it struggles with this case too as it produces different results for float32 and float64.

[mruberry]

Instead of conditioning on CUDA being available, just compare the result to the expected outcome of

tensor([101.,  99., 101.,  99., 101.,  99., 100., 100., 100., 101.])

directly, and add a comment that NumPy produces a different result

[ngimel]

If you change boundary computation to what cpu is doing now (instead of to what cuda is doing) would the results match?

[Aidyn-A]

Yes, the results will exactly match if I slightly change the input tensor or histogram boundaries by a machine epsilon.

[mruberry]

One test tweak needed (see comment inline) and then this should be OK. I'd really like to get the CUDA version of histogramdd added, so we can consistently direct people to that. It's more accurate and more consistent with NumPy.

[kit1980]

@pytorchbot rebase

[kit1980]

@Aidyn-A should we merge this now?

[pytorchmergebot]

@pytorchbot successfully started a rebase job. Check the current status here

[pytorchmergebot]

Successfully rebased fix_consistency_between_cpu_and_cuda_for_histc onto refs/remotes/origin/viable/strict, please pull locally before adding more changes (for example, via git checkout fix_consistency_between_cpu_and_cuda_for_histc && git pull --rebase)

[Aidyn-A]

@kit1980 yes, we can merge it now.
@pytorchbot merge -g

[pytorchmergebot]

Merge started

Your change will be merged once all checks on your PR pass since you used the green (-g) flag (ETA: 0-4 Hours).

Learn more about merging in the wiki.

Questions? Feedback? Please reach out to the PyTorch DevX Team

Advanced Debugging

Check the merge workflow status
here