Counting the points below a coordinate

I have the following problem, for which the solution I use at the moment is too slow.

Instance: a numpy array b of shape (B,2), a sorted numpy array x of shape (X), and a sorted numpy array y of shape (Y).

Note that x = np.unique(b[:,0]) and y = np.unique(b[:,1]), if that makes a difference for the problem.

Task: Build the (X,Y)-array H such that H[i,j] is the number of rows in b whose first entry is less than x[i] and whose second entry is less than y[j].

The following example code solves this:

import numpy as np
b = np.random.random((2000,2))
x = np.unique(b[:,0])
y = np.unique(b[:,1])
H = np.count_nonzero(
    np.logical_and(
        b[:,0,None,None] <= x[None,:,None], 
        b[:,1,None,None] <= y[None,None,:]
    ), 
    axis=0
)

but this gets quite slow if b and thus x and y have a few thousand entries.

How can I do this more efficiently?

source https://stackoverflow.com/questions/77044480/counting-the-points-below-a-coordinate