I am not sure I understand fully your question, but IMO if you have mu1, mu2, sigma_1, sigma_2, why not create two instances of MultivariateNormal using diagonal covariance matrices?
# mu1 [batch_size, n]
# assuming that sigma1 is [batch_size, n] also
cov1 = torch.stack([torch.diag(sigma) for sigma in torch.exp(sigma_1)])
mvn1 = MultivariateNormal(mu1, cov1)

...
# similarly for mvn2

MultivariateNormal will interpret the batch_size as the batch dimension automatically thus mvn1 would have:
batch_shape = batch_size
event_shape = n
sample_shape = ()

when you sample it will take into consideration the batch_shape.
You can compute kl (mvn1, mvn2) using the Pytorch’s implementation.
I hope this can help.