The cross-entropy method will easily allow you to approximate $\mathcal{P}_{q,\epsilon}$ as an ellipsoid, which is likely reasonable if $\epsilon$ is small enough (the$q$ is a global minimum so the hessian is semi definite positive around q$q$)
The idea is to iteratively find a multivariate normal distribution that minimizes its KL-divergence to the distribution $\mathbf{1}_{\mathcal{P}_{q,\epsilon}}$. This will then allow you to efficiently generate random samples from $\mathcal{P}_{q,\epsilon}$.
Note that the C.E method uses KL-divergence, but it has nothing to do with the fact that the problem is about KL-divergence. The answer would be similar for many other types of balls.