算法

您当前的实现很好地进行了向量化，但与可以完成的工作相比，它计算的工作量太大了。主要的问题是这个实现的算法复杂度是**O(n**3)。
这个问题可以用O(n**2)算法解决。这是最佳复杂度**，因为H至少需要填充。关键是对b按y进行排序，以便在计算H行时有效地计算条目的数量。排序后的数组需要通过x进行过滤，这样就不用关心热循环中的x值。注意，这种方法假设x和y是排序的，因为np.unique对值进行排序。

实现

为了提高实现的性能，可以使用Numba和多线程。注意内存布局对算法缓存友好也很重要。
下面是最终的实现：

import numba as nb

@nb.njit('(float64[::1], float64[::1], float64[:,::1])', parallel=True)
def compute(x, y, b):
    H = np.empty((x.size, y.size), dtype=np.int32)

    # Improve the memory layout
    b_x = b[:,0]
    b_y = b[:,1]

    # Sort both b_x and b_y pairs in a way b_y is sorted
    sorted_idx = np.argsort(b_y)
    sorted_b_x = b_x[sorted_idx]
    sorted_b_y = b_y[sorted_idx]

    for i in nb.prange(x.size):
        # Filter by x
        sorted_filtered_b_y = sorted_b_y[sorted_b_x <= x[i]]

        size = sorted_filtered_b_y.size
        count = 0
        idx = 0

        for j in range(y.size):
            cur_y = y[j]
            while idx < size and sorted_filtered_b_y[idx] <= cur_y:
                count += 1
                idx += 1
            H[i, j] = count

    return H

x = np.unique(b[:,0])
y = np.unique(b[:,1])
compute(x, y, b)

性能结果

以下是在我的i5- 9600 KF机器上使用所提供的输入的性能结果：

Naive implementation:    9874.74 ms
This implementation:        3.85 ms

因此，这种实现方式的速度是2565倍！