PS: This article clarifies the relationship between gradient descent and the EM algorithm. Through a unified perspective, it derives ordinary gradient descent, the EM algorithm in pLSA, and the EM algorithm in K-Means. This demonstrates that they are essentially different aspects of the same concept—as the saying goes, “viewed from the front, a ridge; from the side, a peak; different heights and distances yield different views.”
In machine learning, we usually represent the problem we want to solve as a loss function with unknown parameters, such as Mean Squared Error (MSE). We then find ways to solve for the minimum of this function to obtain the optimal parameter values, thereby completing the modeling. Since multiplying a function by -1 turns a maximum into a minimum, we generally refer to finding the minimum. In the field of machine learning, two major directions are commonly taught: 1. Gradient Descent; 2. The EM algorithm (Expectation-Maximization), which is generally used for solving complex maximum likelihood problems.
In typical tutorials, these two methods are described as being completely different, as if they were two competing systems, and the EM algorithm is often described with an air of mystery. However, in reality, both methods are different examples of the same underlying logic. They are “born from the same root” and belong to the same lineage.
Let us begin with the ancient Newton’s method.
Newton’s Iterative Method
Given a complex non-linear function f(x), we wish to find its minimum. Generally, if we assume it is sufficiently smooth, its minimum is also a stationary point where f'(x_0)=0. Thus, the problem is transformed into finding the root of the equation f'(x)=0. For the roots of non-linear equations, we have Newton’s method: x_{n+1} = x_{n} - \frac{f'(x_n)}{f''(x_n)}
However, this approach is detached from geometric meaning and does not allow us to glimpse deeper secrets. We prefer the following logic: at the point x=x_n on the curve y=f(x), we can use an approximate curve to approach the original function. If the approximate curve is easy to minimize, we can use the minimum of this approximate curve to represent the minimum of the original curve:
Clearly, the requirements for the approximate curve are:
1. It must approximate the true curve to some degree; generally, at least first-order approximation is required.
2. It must have a minimum point, and that minimum point must be easy to solve.
Naturally, we can choose a “tangent parabola” for approximation: f(x)\approx g(x) = f(x_n)+f'(x_n)(x-x_n)+\frac{1}{2}f''(x_n)(x-x_n)^2 This parabola has second-order accuracy. For this parabola, the extremum point is: x_n - \frac{f'(x_n)}{f''(x_n)} Thus, we re-derive the iterative formula for Newton’s method: x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)}
If f(x) is sufficiently smooth and has only one global extremum, Newton’s method will converge rapidly (exponentially). However, real-world functions are rarely so ideal, and its disadvantages become apparent:
1. It requires calculating the second derivative, which can be extremely complex for some functions.
2. Since the magnitude of f''(x_n) is uncertain, the opening direction of g(x) is uncertain, so we cannot be sure whether the result is a maximum or a minimum.
Gradient Descent
These two disadvantages are fatal in many problems. Therefore, to solve these issues, we abandon second-order accuracy and replace f''(x_n) with a fixed positive number 1/h: g(x) = f(x_n)+f'(x_n)(x-x_n)+\frac{1}{2h}(x-x_n)^2 This approximate curve has only first-order accuracy, but it eliminates the calculation of the second derivative and ensures that the parabola opens upwards. Therefore, iterating through it guarantees convergence to a minimum (at least a local one). The minimum point of the above g(x) is: x_n - h f'(x_n) So we get the iterative formula: x_{n+1} = x_n - h f'(x_n) For high-dimensional spaces, this is: \boldsymbol{x}_{n+1} = \boldsymbol{x}_n - h \nabla(\boldsymbol{x}_n) This is the famous Gradient Descent method. Of course, it has its own problems, but many improved algorithms, such as Stochastic Gradient Descent, revolve around it.
Here, we understand Gradient Descent as the result of using an approximate parabola. Seeing it this way, readers might wonder: Why must I use a parabola? Can’t I use other curves? Of course you can. For many problems, Gradient Descent might even complicate things—meaning the parabola fails. In such cases, we must consider other forms of approximation. In fact, other approximation schemes are basically referred to as the EM algorithm, which happens to exclude only its sibling, Gradient Descent, which is quite puzzling.
Maximum Likelihood
When estimating probabilities, our optimization goal is usually the maximum likelihood function rather than MSE. For example, if we are building a language model, we need to estimate the co-occurrence probability p(x,y) of any two words x,y. Suppose in a corpus of size N, x,y co-occur \#(x,y) times; we then get the statistical result: \tilde{p}(x,y)=\frac{\#(x,y)}{N} This is a basic result, but it suffers from sparsity issues, and storing a result for every possible word pair would consume too much memory.
A better solution is to assume p(x,y) can be represented by a function p(x,y;\theta) with unknown parameters \theta (where \theta might be a vector), such as a neural network, and then find a way to optimize the parameters. So, what is the optimization goal? Note that if we use MSE, the biggest problem is that we cannot guarantee the result will be non-negative, whereas probabilities must be non-negative.
For probability problems, statisticians proposed a more natural scheme—the maximum likelihood function. Philosophers often say “existence is reasonable”; the idea of the maximum likelihood function is even more thorough: “existence is the most reasonable.” If x,y co-occur \#(x,y) times, since this event occurred, it must be the most reasonable. Therefore, the probability function: \prod_{x,y} p(x,y;\theta)^{\#(x,y)} should be maximized. Taking the logarithm gives: \sum_{x,y} \#(x,y)\log p(x,y;\theta) We should maximize this function, which is called the log-likelihood function. Clearly, \#(x,y) can be replaced by the statistical frequency \tilde{p}(x,y), and the result is equivalent: \sum_{x,y} \tilde{p}(x,y)\log p(x,y;\theta) In fact, multiplying it by -1 gives: S=-\sum_{x,y} \tilde{p}(x,y)\log p(x,y;\theta) We give this a special name—Cross-Entropy, which is a common loss function in machine learning. That is, maximizing the likelihood function is equivalent to minimizing the cross-entropy. If we do not pre-specify the form of p(x,y;\theta) but estimate p(x,y;\theta) directly, it is not hard to calculate that p(x,y;\theta)=\tilde{p}(x,y), which is exactly what we expect and demonstrates the rationality of the maximum likelihood function as an optimization goal.
EM Algorithm
For the optimization of cross-entropy, we usually try Gradient Descent. However, in many cases, Gradient Descent is not effective, and we prefer to use the EM algorithm, which is called the “Algorithm of God.”
Using the previous language model as an example, to improve generalization, we transform p(x,y) into p(x|y)p(y) and decompose p(x|y) as p(x|y)=\sum_z p(x|z)p(z|y). The significance of this decomposition was mentioned in the article SVD Decomposition (II): Why SVD Implies Clustering?. Here, z can be understood as a category or topic; p(x|y) is the probability of x following y, p(z|y) is the probability that y belongs to topic z, and p(x|z) is the probability of x appearing in topic z. Generally, the number of z is much smaller than the number of x,y, thereby reducing the total number of parameters.
At this point, the cross-entropy is: S=-\sum_{x,y} \tilde{p}(x,y)\log \sum_z p(x|z)p(z|y)p(y) We can assume p(y) can be accurately estimated by \tilde{p}(y), so this term only contributes a constant. Thus, the equivalent optimization goal is: S=-\sum_{x,y} \tilde{p}(x,y)\log \sum_z p(x|z)p(z|y) Here p(x|z) and p(z|y) are the parameters to be solved (traversing all possible x,y,z combinations).
Its gradients are: \begin{aligned}&\frac{\partial S}{\partial p(x|z)} = -\sum_{y} \frac{\tilde{p}(x,y)}{\sum_z p(x|z)p(z|y)}p(z|y)\\ &\frac{\partial S}{\partial p(z|y)} = -\sum_{x} \frac{\tilde{p}(x,y)}{\sum_z p(x|z)p(z|y)}p(x|z)\end{aligned}
Direct Gradient Descent is unfeasible because p(x|z) and p(z|y) must be non-negative and satisfy the constraints: \sum_x p(x|z) = 1,\quad \sum_z p(z|y)=1 Gradient Descent cannot guarantee non-negativity constraints, which dooms it to be ineffective for such problems. Recalling the logic of deriving Gradient Descent at the beginning of the article—using an approximate curve to iterate—Gradient Descent uses a parabolic approximation. Here, we do not use a parabolic approximation (as it cannot guarantee positive definiteness), but instead construct a new approximation. Suppose we have performed n iterations and obtained estimates p_n(x|z) and p_n(z|y). According to the gradient formula, the current gradients are: \begin{aligned}&\frac{\partial S}{\partial p_n(x|z)} = -\sum_{y} \frac{\tilde{p}(x,y)}{\sum_z p_n(x|z)p_n(z|y)}p_n(z|y)\\ &\frac{\partial S}{\partial p_n(z|y)} = -\sum_{x} \frac{\tilde{p}(x,y)}{\sum_z p_n(x|z)p_n(z|y)}p_n(x|z)\end{aligned}
The difficulty of the original problem lies in the summation inside the \log. If we could move the summation outside, it would be simple. Therefore, we might consider an approximate function like this: S_n'=-\sum_{x,y} \tilde{p}(x,y)\sum_z C_{x,y,z,n} \log [p(x|z)p(z|y)] where C is a constant (in the current iteration). In this way, S_n' also has a minimum value that can be solved exactly. Naturally, we want the gradient of S' to be the same as the gradient of the original S (possessing first-order accuracy). The gradient of S' is: \begin{aligned}&\frac{\partial S'}{\partial p_n(x|z)} = -\sum_{y} \frac{\tilde{p}(x,y)C_{x,y,z,n}}{p_n(x|z)}\\ &\frac{\partial S'}{\partial p_n(z|y)} = -\sum_{x} \frac{\tilde{p}(x,y)C_{x,y,z,n}}{p_n(z|y)}\end{aligned} Comparing the two sets of gradients, we get: C_{x,y,z,n}=\frac{p_n(x|z)p_n(z|y)}{\sum_z p_n(x|z)p_n(z|y)} That is to say, once we have a set of initial parameters, we can substitute them into the above equation to find C_{x,y,z,n}, and then find the parameters that minimize S_n' to serve as p_{n+1}(x|z) and p_{n+1}(z|y), and iterate accordingly.
This part basically covers the solution process of the pLSA model. For more details, one can refer to tutorials on Probabilistic Latent Semantic Analysis. For this article, the point is to show that the EM algorithm, like Gradient Descent, comes from the same source—both are based on approximate curve fitting. It is not some mystical method, and the search for this approximate function is well-founded. Almost all online tutorials directly give the expression for S' (usually called the Q-function) in a way that feels like metaphysics, which I find quite off-putting.
K-Means
K-Means clustering is easy to understand: given the coordinates of N points \boldsymbol{x}_i,\,i=1,\dots,N, we find a way to divide these points into K categories. Each category has a cluster center \boldsymbol{c}_j,\,j=1,\dots,K. Naturally, the category a point belongs to is the one represented by the nearest cluster center \boldsymbol{c}_j, where distance is defined as Euclidean distance.
Therefore, the main task of K-Means clustering is to find the cluster centers \boldsymbol{c}_j. We naturally hope that each cluster center is exactly at the “center” of the category. Expressed as a function, we want to minimize the following function L: L=\sum_{i=1}^N \min\bigg\{|\boldsymbol{x}_i-\boldsymbol{c}_1|^2,|\boldsymbol{x}_i-\boldsymbol{c}_2|^2,\dots,|\boldsymbol{x}_i-\boldsymbol{c}_K|^2\bigg\} where the \min operation ensures that each point belongs only to the class closest to it.
If we directly use Gradient Descent to optimize L, we will encounter great difficulties. This is not because the \min operation is hard to differentiate, but because this is an NP problem where the theoretical convergence time grows exponentially with N. In this case, we also use the EM algorithm, which manifests as:
1. Randomly select K points as initial cluster centers;
2. Given the K cluster centers, determine which category each point belongs to, and then use the average coordinates of all points in the same category as the new cluster center.
This method usually converges in a few iterations. So, what is the reason for doing this?
We still follow the logic of approximate curve fitting. But the problem is, what if the current \min is not differentiable? We can consider a smooth approximation and then take the limit. The answer is here: Seeking a Smooth Maximum Function. Taking a sufficiently large M, we can assume (the minimum equals the negative of the maximum of the negatives): \begin{aligned}&\min\bigg\{|\boldsymbol{x}_i-\boldsymbol{c}_1|^2,|\boldsymbol{x}_i-\boldsymbol{c}_2|^2,\dots,|\boldsymbol{x}_i-\boldsymbol{c}_K|^2\bigg\}\\ =&-\frac{1}{M}\ln\bigg(e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_1|^2}+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_2|^2}+\dots+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_K|^2}\bigg)\end{aligned} Then: L=-\sum_{i=1}^N \frac{1}{M}\ln\bigg(e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_1|^2}+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_2|^2}+\dots+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_K|^2}\bigg) Now we can find the gradient: \frac{\partial L}{\partial \boldsymbol{c}_j}=\sum_{i=1}^N \frac{2e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_j|^2 } }{e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_1|^2}+\dots+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_K|^2}}(\boldsymbol{c}_j-\boldsymbol{x}_i) Let the result of the n-th iteration be \boldsymbol{c}^{(n)}_j; the gradient for this round is: \frac{\partial L}{\partial \boldsymbol{c}^{(n)}_j}=\sum_{i=1}^N \frac{2e^{-M|\boldsymbol{x}_i-\boldsymbol{c}^{(n)}_j|^2 } }{e^{-M|\boldsymbol{x}_i-\boldsymbol{c}^{(n)}_1|^2}+\dots+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_K^{(n)}|^2}}(\boldsymbol{c}^{(n)}_j-\boldsymbol{x}_i) Based on the characteristics of this formula, we can look for an approximate curve (or rather, a hypersurface): L'=\sum_{i=1}^N \sum_{j=1}^K C^{(n)}_{i,j} |\boldsymbol{x}_i-\boldsymbol{c}_j|^2 where C^{(n)}_{i,j} is to be determined and is a constant in each round. This is just a quadratic function, and it is easy to find that its minimum point is: \boldsymbol{c}_j = \frac{\sum_{i=1}^N C^{(n)}_{i,j}\boldsymbol{x}_i}{\sum_{i=1}^N C^{(n)}_{i,j}} which is the weighted average of \boldsymbol{x}_i.
Similarly, we want this approximate curve to have at least first-order accuracy compared to the original function, so we find its derivative: \frac{\partial L'}{\partial \boldsymbol{c}_j}=\sum_{i=1}^N 2C^{(n)}_{i,j} (\boldsymbol{c}_j-\boldsymbol{x}_i) Comparing it with the derivative of the original function, we easily get: C^{(n)}_{i,j} = \frac{e^{-M|\boldsymbol{x}_i-\boldsymbol{c}^{(n)}_j|^2 } }{e^{-M|\boldsymbol{x}_i-\boldsymbol{c}^{(n)}_1|^2}+\dots+e^{-M|\boldsymbol{x}_i-\boldsymbol{c}_K^{(n)}|^2}} At this point, we get the iterative formula: \boldsymbol{c}^{(n+1)}_j = \frac{\sum_{i=1}^N C^{(n)}_{i,j}\boldsymbol{x}_i}{\sum_{i=1}^N C^{(n)}_{i,j}}
So far, we have derived all the steps, but we are using a continuous approximation. Finally, we take the limit as M\to\infty. After taking the limit, the problem simplifies. From the above formula, we can derive: \lim_{M\to\infty} C^{(n)}_{i,j} = \Delta^{(n)}_{i,j} = \left\{\begin{aligned}&1, \text{ if for a fixed } i, \text{ the distance from } j \text{ to } i \text{ is minimal}\\ &0, \text{ otherwise}\end{aligned}\right. Simply put, if we view \Delta^{(n)}_{i,j} as an N \times K matrix, each row of the matrix can have only one 1, and the rest are 0. If the j-th element in the i-th row is 1, it means the cluster center closest to \boldsymbol{x}_i is \boldsymbol{c}_j. At this point, the iterative formula becomes: \boldsymbol{c}^{(n+1)}_j = \frac{\sum_{i=1}^N \Delta^{(n)}_{i,j}\boldsymbol{x}_i}{\sum_{i=1}^N \Delta^{(n)}_{i,j}} According to the meaning of \Delta^{(n)}_{i,j}, this is nothing more than saying:
\boldsymbol{c}^{(n+1)}_j is the average of those points closest to \boldsymbol{c}^{(n)}_j.
This yields the iterative algorithm commonly used to solve the K-Means problem, which is also known as the EM algorithm.
Summary
As we can see, the so-called EM algorithm is not a specific method, but a class of methods, or a strategy. Gradient Descent is originally a special case within this class of methods; they are entirely the same thing, and strictly speaking, it should not be excluded. The so-called “Algorithm of God” is actually just an iterative method. Through its own iterative updates, it can approach perfection (the optimal solution), much like biological evolution. It is as if it were carefully designed by a creator, so it is no wonder it is called the Algorithm of God.
When reposting, please include the original address: https://kexue.fm/archives/4277
For more detailed reposting matters, please refer to: Scientific Space FAQ