How to use KL divergence to compare two distributions?

I am trying to model the probability distribution of a multi-dimensional dataset where all the values are discrete. Suppose the training data (represented by T) is of the shape (m, n) where n is the number of features and m is the number of samples. I use a neural network to generate a dataset G which is also of shape (m, n).

I calculate the empirical $$pmfs$$ for each of the features separately: $$P = (p_{i1}, p_{i2}, ... p_{in})$$ representing the actual distribution and $$Q = (q_{i1}, q_{i2}, ..., q_{in})$$ representing the generated distribution. Here each $$p_{i}$$ and $$q_{i}$$ is multi-dimensional since each feature can take multiple categorical values. Finally, I am calculating the KLD as:

$$KLD = \frac{1}{n} \times q_{i}log(\frac{q_{i}}{p_{i}})$$

The logs and multiplication are done element-wise, since $$p_{i}$$ and $$q_{i}$$ are multi-dimensional.

While training, the objective is to minimize the above sum. Is the above calculation correct?

Another method would be to calculate the $$P$$ and $$Q$$ for the entire distribution together, but that seems like a bad idea since each feature can take a number of values, and the training sample size is pretty small (around 5,000 samples).

How do I calculate KL-divergence between two multidimensional distributions? is a somewhat similar question but there was no clear answer on that thread.