Your kl is not a KL divergence between q and p.
If you want to estimate KL divergence between q and p using sampling, then you just need to get samples from q and then calculate log q(sample) - log p(samples) and take average over number of samples. But computing KL between two relaxed onehots suffers from numerical instability, it can be better to compute KL between ExpRelaxedOneHotCategorical as KL is invariant under invertible transformations. For more information you can read https://arxiv.org/pdf/1611.00712.pdf 34