In [1]: import numpy as np
In [2]: a = np.random.rand(25*10**6)
In [3]: b = np.random.rand(2*10**5)
In [4]: %timeit np.isin(a,b)
3.88 s ± 119 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

like so?

Also, I was asking what do you do after you get the mask – it may not be the best thing to do depends on what your ultimate goal is.

In [1]: import numpy as np
In [2]: a = np.random.rand(25*10**6)
In [3]: b = np.random.rand(2*10**5)
In [4]: %timeit np.isin(a,b)
3.88 s ± 119 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

edit: they sort it:

in Julia you can also just use insorted(), notice the time includes the time to sort b

julia> @be insorted.(a, Ref(sort(b))) evals=10
Benchmark: 1 sample with 10 evaluations
2.566 s (13 allocs: 6.040 MiB, 0.02% gc time)