The KL Divergence is quite easy to compute in closed form for simple distributions -such as Gaussians- but has some not-very-nice properties. For example, it is not symmetrical (thus it is not a metric) and it does not respect the triangular inequality.
What is the reason it is used so often in ML? Aren't there other statistical distances that can be used instead?
In ML we always deal with unknown probability distributions from which the data comes. The most common way to calculate the distance between real and model distribution is $KL$ divergence.
Although there are other loss functions (e.g. MSE, MAE), $KL$ divergence is natural when we are dealing with probability distributions. It is a fundamental equation in information theory that quantifies, in bits, how close two probability distributions are. It is also called relative entropy and, as the name suggests, it is closely related to entropy, which in turn is a central concept in information theory. Let's recall the definition of entropy for a discrete case:
$$ H = -\sum_{i=1}^{N} p(x_i) \cdot \text{log }p(x_i) $$
As you observed, entropy on its own is just a measure of a single probability distribution. If we slightly modify this formula by adding a second distribution, we get $KL$ divergence:
$$ D_{KL}(p||q) = \sum_{i=1}^{N} p(x_i)\cdot (\text{log }p(x_i) - \text{log }q(x_i)) $$
where $p$ is a data distribution and $q$ is model distribution.
As we can see, $KL$ divergence is the most natural way to compare 2 distributions. Moreover, it's pretty easy to calculate. This article provides more intuition on this:
Essentially, what we're looking at with the KL divergence is the expectation of the log difference between the probability of data in the original distribution with the approximating distribution. Again, if we think in terms of $log_2$ we can interpret this as "how many bits of information we expect to lose".
Cross-entropy is commonly used in machine learning as a loss function where we have softmax (or sigmoid) output layer, since it represents a predictive distribution over classes. The one-hot output represents a model distribution $q$, while true labels represent a target distribution $p$. Our goal is to push $q$ to $p$ as close as possible. We could take a mean squared error over all values, or we could sum the absolute differences, but the one measure which is motivated by information theory is cross-entropy. It gives the average number of bits needed to encode samples distributed as $p$, using $q$ as the encoding distribution.
Cross-entropy based on entropy and generally calculates the difference between two probability distributions and closely related to $KL$ divergence. The difference is that it calculates the the total entropy between the distributions, while $KL$ divergence represents relative entropy. Corss-entropy can be defined as follows:
$$ H(p, q) = H(p) + D_{KL}(p \parallel q) $$
The first term in this equation is the entropy of the true probability distribution $p$ that is omitted during optimization, since the entropy of $p$ is constant. Hence, minimizing cross-entropy is the same as optimizing $KL$ divergence.
It can be also shown that maximizing the (log) likelihood is equivalent to minimizing the cross entropy.
As you mentioned, $KL$ divergence is not symmetrical. But in most cases this is not critical, since we want to estimate the model distribution by pushing it towards real one, but not vice versa. There is also a symmetrized version called Jensen–Shannon divergence: $$ D_{JS}(p||q)=\frac{1}{2}D_{KL}(p||m)+\frac{1}{2}D_{KL}(q||m) $$ where $m=\frac{1}{2}(p+q)$.
The main disadvantage of $KL$ is that both the unknown distribution and the model distribution must have support. Otherwise the $D_{KL}(p||q)$ becomes $+\infty$ and $D_{JS}(p||q)$ becomes $log2$
Second, it should be noted that $KL$ is not a metric, since it violates triangle inequality. That is, in some cases it won't tell us if we are going the right direction when estimating our model distribution. Here is an example taken from this answer. Given two discrete distributions $p$ and $q$, we calculate $KL$ divergence and Wasserstein metric:
As you can see, $KL$ divergence remained the same, while the Wasserstein metric decreased.
But as mentioned in comments, Wasserstein metric is highly intractable in a continuous space. We still can use it by applying the Kantorovich-Rubinstein duality used in Wasserstein GAN. You can also find more on this topic in this article.
The 2 drawbacks of $KL$ can be mitigated by adding noise. More on it in this paper
This question is very general in the sense that the reason may differ depending on the area of ML you are considering. Below are two different areas of ML where the KL-divergence is a natural consequence:
Let $X_1, \dots, X_n$ be i.i.d. Bernoulli RVs with parameter $p$. $$P(\sum_i X_i \geq (p+t)n) \leq \inf_\lambda M_X (\lambda) e^{-\lambda t} = \exp(-n D_{KL}(p+t||p)).$$
