Remove pixel patch in image which is stored as array
On my computer it takes 1.745 seconds to run the code in the post.
There's no need for the array of random indexes to be two-dimensional:

IDX = np.random.randint(0,P,(N,1))


In fact this is harmful for performance, because it means that x0[i] is an array of length 1 (not a scalar) and so img[x0[i]+x,y0[i]+y] requires "fancy indexing" which is slower than normal indexing.

It would be simpler to make the array of indexes one-dimensional:

IDX = np.random.randint(P, size=N)


This reduces the runtime to about 0.459 seconds (26.3% of the original).
There is no need to reassign I[i] at the end of the loop. When you call the reshape method on a NumPy array, what you get is a view onto the original array (not a copy) if possible. (And it is possible in this case.) So updating the view also updates the original.

This reduces the runtime to about 0.449 seconds (25.8%).
Instead of looping over range(N) and then looking up I[i] and x0[i] and y0[i], use zip to loop over all the arrays simultaneously:

for img, xx, yy in zip(I, x0, y0):
    img = img.reshape(q,q)
    for x in range(ps):
        for y in range(ps):
            if xx + x &lt; q and yy + y &lt; q:
                img[xy + x, yy + y] = 0.0


This reduces the runtime to about 0.358 seconds (20.5%).
Instead of looping over all the pixels in the patch and updating each pixel individually, use slices to update the whole region in one step:

for image, x, y in zip(I, x0, y0):
    image.reshape(q, q)[x:x + ps, y:y + ps] = 0.0


This works because NumPy (and Python generally) ensures that the bounds of a slice do not go beyond the end of the array. See the slicing documentation:


  The slice of \$s\$ from \$i\$ to \$j\$ is defined as the sequence of items with index \$k\$ such that \$i \le k &lt; j\$. If \$i\$ or \$j\$ is greater than len(s), use len(s).


This reduces the runtime to about 0.025 seconds (1.4%).
We can vectorize the additions x + ps and y + ps:

for image, x, y, x1, y1 in zip(I, x0, y0, x0 + ps, y0 + ps):
    image.reshape(q, q)[x:x1, y:y1] = 0.0


This reduces the runtime to about 0.021 seconds (1.2%).
We could avoid the reshape inside the loop by doing a single reshape of the whole I array:

images = I.reshape(N, q, q)


and then:

for image, x, y, x1, y1 in zip(images, x0, y0, x0 + ps, y0 + ps):
    image[x:x1, y:y1] = 0.0


This reduces the runtime to about 0.018 seconds (1.0%).
We can halve the number of indexing operations by indexing the images array just once on each loop iteration:

for i, x, y, x1, y1 in zip(range(N), x0, y0, x0 + ps, y0 + ps):
    images[i, x:x1, y:y1] = 0.0


This reduces the runtime to about 0.011 seconds (0.6%).


That's about 150 times speedup overall, so calling this a million times will still take about 3 hours on my computer. There may be other improvements to be had if only we could see more of your code, but you'll need to make a new post for that.