Possible bug in the implementation of histogram in histogram based gradient boosting.

### Describe the bug

I was looking at the code in the file `sklearn/ensemble/_hist_gradient_boosting/histogram.pyx` and I believe there might be a bug. The note at the top of the file describes the following:

```
# Notes:
# - IN views are read-only, OUT views are write-only
# - In a lot of functions here, we pass feature_idx and the whole 2d
#   histograms arrays instead of just histograms[feature_idx]. This is because
#   Cython generated C code will have strange Python interactions (likely
#   related to the GIL release and the custom histogram dtype) when using 1d
#   histogram arrays that come from 2d arrays.
# - The for loops are un-wrapped, for example:
#
#   for i in range(n):
#       array[i] = i
#
#   will become
#
#   for i in range(n // 4):
#       array[i] = i
#       array[i + 1] = i + 1
#       array[i + 2] = i + 2
#       array[i + 3] = i + 3
#
#   This is to hint gcc that it can auto-vectorize these 4 operations and
#   perform them all at once.
```

This is used in several places in the file for example here:

```

cpdef void _build_histogram(
        const int feature_idx,
        const unsigned int [::1] sample_indices,  # IN
        const X_BINNED_DTYPE_C [::1] binned_feature,  # IN
        const G_H_DTYPE_C [::1] ordered_gradients,  # IN
        const G_H_DTYPE_C [::1] ordered_hessians,  # IN
        hist_struct [:, ::1] out) nogil:  # OUT
    """Return histogram for a given feature."""
    cdef:
        unsigned int i = 0
        unsigned int n_node_samples = sample_indices.shape[0]
        unsigned int unrolled_upper = (n_node_samples // 4) * 4

        unsigned int bin_0
        unsigned int bin_1
        unsigned int bin_2
        unsigned int bin_3
        unsigned int bin_idx

    for i in range(0, unrolled_upper, 4):
        bin_0 = binned_feature[sample_indices[i]]
        bin_1 = binned_feature[sample_indices[i + 1]]
        bin_2 = binned_feature[sample_indices[i + 2]]
        bin_3 = binned_feature[sample_indices[i + 3]]

        out[feature_idx, bin_0].sum_gradients += ordered_gradients[i]
        out[feature_idx, bin_1].sum_gradients += ordered_gradients[i + 1]
        out[feature_idx, bin_2].sum_gradients += ordered_gradients[i + 2]
        out[feature_idx, bin_3].sum_gradients += ordered_gradients[i + 3]

        out[feature_idx, bin_0].sum_hessians += ordered_hessians[i]
        out[feature_idx, bin_1].sum_hessians += ordered_hessians[i + 1]
        out[feature_idx, bin_2].sum_hessians += ordered_hessians[i + 2]
        out[feature_idx, bin_3].sum_hessians += ordered_hessians[i + 3]

        out[feature_idx, bin_0].count += 1
        out[feature_idx, bin_1].count += 1
        out[feature_idx, bin_2].count += 1
        out[feature_idx, bin_3].count += 1

    for i in range(unrolled_upper, n_node_samples):
        bin_idx = binned_feature[sample_indices[i]]
        out[feature_idx, bin_idx].sum_gradients += ordered_gradients[i]
        out[feature_idx, bin_idx].sum_hessians += ordered_hessians[i]
        out[feature_idx, bin_idx].count += 1
```

However, there is an issue. Here we are indexing over bins for example `bin_0`, `bin_1`, and so on. Nevertheless, it is possible that the bins are the same (this is in fact not too unlikely because there are only 256 possible bins). The issue is that if this happens, and it is the case that the compiler executes the function as the note describes, then we will get race conditions because it is possible that we modify the same memory in `out` simultaneously because of the incorrect parallel execution.  I am not sure if this is a bug or not because it depends on how this is compiled -- but at the very least the code and the comments are incorrect. 

### Steps/Code to Reproduce

There is no code to reproduce the problem. 

### Expected Results

The expected result is that the computation should have the same output as done without vectorization and roop unrolling but I am not sure this is the case. 

### Actual Results

It is unclear what the actual results are but there is an issue at the very least with what is written in the code. 

### Versions

It is present in branch 1.2.x of the repo (and the main branch as of Dec 22 2022)

COMMENT:
@velezbeltran It is astonishing and very rare to get a bug report by pure looking at the code! May I ask how you came to look that deep into the code?

To your question, I agree, the comment could be improved. Concerning the write collisions, I can give the all-clear. The code itself is fine and has no bugs or problems with access patterns as it is to be read sequentially. The compiler is only allowed to do code transformations that give the same output.
The point is, there is hardly any vectorization happening. See the analysis of the more or less equivalent C code here: https://godbolt.org/z/af6671xMa.

Lately, I played with calculating histograms, see https://lorentzen.ch/index.php/2022/10/31/histograms-gradient-boosted-trees-group-by-queries-and-one-hot-encoding/. This issue and the findings from godbolt made me think if there are opportunities for improvement and there might be some. So thanks for this issue😄

Labels: Bug, Needs Triage