【Searching for Text】· (2) From MCMC to Simulated Annealing
In the previous article, we introduced the concept of "constrained text generation," pointing out that certain conditional text generation tasks can be accomplished in an unsupervised way by quantifying an objective and sampling from it. That article also introduced two methods, "importance sampling" and "rejection sampling," and pointed out that for high-dimensional spaces, the easy-to-sample distributions that these methods rely on are often hard to design, making it difficult for them to satisfy our sampling needs.
At this point, we need to bring in one of the most important algorithms in the sampling world, "Markov Chain Monte Carlo (MCMC)." It combines Markov chains with the Monte Carlo method, making it possible (at least in theory) to sample from many high-dimensional distributions. It is also one of the essential foundational algorithms for the constrained text generation applications we will introduce later. This article attempts to give a basic introduction to it.
Markov Chains
A Markov chain is essentially a "memoryless" random walk process. Based on a transition probability $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$, starting from an initial state $\boldsymbol{x}_0$, at each step the next state is chosen randomly according to this transition probability, thereby forming a random sequence of states $\boldsymbol{x}_0, \boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_t, \cdots $. We want to examine, for a sufficiently large number of steps $t$, the distribution followed by $\boldsymbol{x}_t$ — that is, the "stationary distribution" of this Markov chain.
$\setCounter{11}$Suppose that the distribution followed by $\boldsymbol{x}_t$ is $p^{(t)}(\boldsymbol{x}_t)$. Then the distribution followed by $\boldsymbol{x}_{t+1}$ is
\begin{equation}p^{(t+1)}(\boldsymbol{x}_{t+1}) = \sum_{\boldsymbol{x}_t} p(\boldsymbol{x}_{t+1} \leftarrow \boldsymbol{x}_t)p^{(t)}(\boldsymbol{x}_t)\end{equation}
If equilibrium has already been reached, then we should have $p^{(t+1)}=p^{(t)}$. So if $p(\boldsymbol{x})$ is the stationary distribution, it satisfies
\begin{equation}p(\boldsymbol{y}) = \sum_{\boldsymbol{x}} p(\boldsymbol{y} \leftarrow \boldsymbol{x})p(\boldsymbol{x})\label{eq:stable}\end{equation}
In other words, the stationary distribution is a nonzero solution of the above equation with respect to $p(\boldsymbol{x})$.
A question that then needs answering is: is the nonzero solution of the above equation unique? The answer to this question determines whether the stationary distribution of a Markov chain depends on the initial state $\boldsymbol{x}_0$. Unfortunately, the answer is: not necessarily. Briefly speaking, the condition for the nonzero solution to be unique is that any two states of the Markov chain are connected — more concretely, for any two states $\boldsymbol{x},\boldsymbol{y}$, there exists a chain of states $\boldsymbol{x}=\boldsymbol{x}_0,\boldsymbol{x}_1,\cdots,\boldsymbol{x}_{n-1},\boldsymbol{x}_n=\boldsymbol{y}$ such that
\begin{equation}\begin{aligned} &p(\boldsymbol{x}_{1} \leftarrow \boldsymbol{x}_0) > 0 \\ &p(\boldsymbol{x}_{2} \leftarrow \boldsymbol{x}_1) > 0 \\ &\qquad\vdots \\ &p(\boldsymbol{x}_{n} \leftarrow \boldsymbol{x}_{n-1}) > 0 \end{aligned}\end{equation}
Put plainly, this means that the probability of jumping from state $\boldsymbol{x}$ to $\boldsymbol{y}$ (not necessarily in one step — multiple steps are allowed) is greater than 0. This criterion is actually not hard to understand intuitively: it means that the condition for a unique stationary distribution is that the Markov chain must not contain any "isolated islands." Otherwise, if $\boldsymbol{x}_0$ is inside an island, it will wander forever within that island, and if $\boldsymbol{x}_0$ is outside the island, it will wander forever outside it — so different initial states would lead to different stationary distributions.
Detailed Balance
Since a Markov chain has a stationary distribution, and constructing a Markov chain only requires a transition matrix, if we could construct an easy-to-sample transition matrix whose stationary distribution happens to be the distribution $p(\boldsymbol{x})$ we want to sample from, then after enough iterations, wouldn't that be equivalent to sampling from $p(\boldsymbol{x})$? This is the idea behind all MCMC methods. Of course, a Markov chain often needs many iterations before reaching stationarity, so sampling tends to be time-consuming; but even so, it can eventually complete the sampling process within a finite cost, and so remains practical in many situations.
Assuming MCMC is feasible, the key question is clearly how to construct a transition probability $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$ whose stationary distribution is our given distribution $p(\boldsymbol{x})$. At the same time, as explained in the previous article, we typically only know a quantity $\rho(\boldsymbol{x})$ proportional to $p(\boldsymbol{x})$, without knowing the normalization factor. Therefore, when constructing $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$, we cannot rely on the normalization factor of $p(\boldsymbol{x})$.
For this purpose, we need to use a "detailed balance condition":
Detailed balance condition: If the distribution $p(\boldsymbol{x})$ and the transition probability $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$ satisfy the identity \begin{equation}p(\boldsymbol{y}\leftarrow\boldsymbol{x})p(\boldsymbol{x})=p(\boldsymbol{x}\leftarrow\boldsymbol{y})p(\boldsymbol{y})\end{equation},
then $p(\boldsymbol{x})$ is the stationary distribution of $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$.
Although this sounds sophisticated, it is actually very easy to prove: just sum both sides over $\boldsymbol{x}$ and we immediately obtain equation $\eqref{eq:stable}$; satisfying equation $\eqref{eq:stable}$ means that $p(\boldsymbol{x})$ is the stationary distribution of $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$. The "detailed balance condition" is a sufficient but not necessary condition for identifying a stationary distribution; its value lies in providing a convenient way to construct transition probabilities for an arbitrary distribution.
Transition Probability
Given the stationary distribution $p(\boldsymbol{x})$ and an arbitrary reference transition probability $q(\boldsymbol{y}\leftarrow\boldsymbol{x})$, if it satisfies the detailed balance condition, then all is well — we can just use the reference transition probability $q(\boldsymbol{y}\leftarrow\boldsymbol{x})$ directly as the final probability. If it doesn't, let us denote $\alpha(\boldsymbol{x}\leftarrow\boldsymbol{y})=q(\boldsymbol{y}\leftarrow\boldsymbol{x})p(\boldsymbol{x})$; then we have the following identity, which obviously holds:
\begin{equation}\underbrace{\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})q(\boldsymbol{y}\leftarrow\boldsymbol{x})}_{\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})}p(\boldsymbol{x})=\underbrace{\alpha(\boldsymbol{x}\leftarrow\boldsymbol{y})q(\boldsymbol{x}\leftarrow\boldsymbol{y})}_{\tilde{q}(\boldsymbol{x}\leftarrow\boldsymbol{y})}p(\boldsymbol{y})\label{eq:feihua}\end{equation}
In fact, this identity is nearly a tautology — it's equivalent to saying "although $a\neq b$, $ab=ba$ still always holds." But it is precisely this seemingly trivial identity that gives us the insight: if we use $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ as the transition probability, wouldn't the corresponding stationary distribution be exactly $p(\boldsymbol{x})$?
It's not quite that simple, however, because the $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ constructed this way is generally not normalized, and being unnormalized means it isn't itself a valid transition probability, in which case there's no such thing as a stationary distribution to speak of. Note that although we may not need to know the normalization factor in order to actually perform sampling, when it comes to the theoretical distribution, we do need normalization to hold. Some readers might wonder: can't we just manually divide by a normalization factor? No, because after dividing by the normalization factor, the detailed balance condition may no longer hold — even equation $\eqref{eq:stable}$ might fail. The author looked through many reference materials on this topic, none of which clarified this issue; they all simply used $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ directly as the transition probability, which left me puzzled for quite a while.
In fact, MCMC methods don't normalize $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ at all — instead, they shift all of the "leftover" probability onto the state itself. In mathematical terms, the transition probability ultimately used is actually:
\begin{equation}p(\boldsymbol{y}\leftarrow\boldsymbol{x})=\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x}) + \left(1 - \sum_{\boldsymbol{y}} \tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})\right)\delta(\boldsymbol{y}\leftarrow\boldsymbol{x}) \end{equation}
where $\delta(\boldsymbol{y}\leftarrow\boldsymbol{x}) = \left\{\begin{aligned}1,\boldsymbol{y}=\boldsymbol{x} \\ 0, \boldsymbol{y}\neq \boldsymbol{x}\end{aligned}\right.$ represents the transition probability of staying in the same state (i.e., "never changing"). So the definition of $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$ is really just piling all the excess probability onto the "no change" operation, which means it obviously satisfies normalization. As for the detailed balance condition, it can be verified by substitution, and it also holds — largely because for any $p(\boldsymbol{x})$ and $f(\boldsymbol{x},\boldsymbol{y})$, the following identity always holds:
\begin{equation}f(\boldsymbol{x},\boldsymbol{y})\delta(\boldsymbol{y}\leftarrow \boldsymbol{x})p(\boldsymbol{x})=f(\boldsymbol{y},\boldsymbol{x})\delta(\boldsymbol{x}\leftarrow \boldsymbol{y})p(\boldsymbol{y})\end{equation}
So $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$ is the true transition probability we are actually seeking. If the state space is finite, then the transition probability corresponds to a finite matrix, and the above result essentially says: "normalization can be achieved by adjusting the diagonal elements of the transition matrix, without affecting whether the detailed balance condition holds."
The MCMC Method
How do we implement sampling from the above $p(\boldsymbol{y}\leftarrow\boldsymbol{x})$? It's simple: first implement sampling from $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$, and then, if the condition for sampling from $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ is not satisfied, just stay put; and $\tilde{q}(\boldsymbol{y}\leftarrow\boldsymbol{x})=\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})q(\boldsymbol{y}\leftarrow\boldsymbol{x})$, we assume that sampling from $q(\boldsymbol{y}\leftarrow\boldsymbol{x})$ is easy. Recalling the rejection sampling introduced in the previous article, we can further sample a $\varepsilon\sim U[0,1]$, and via $ \alpha(\boldsymbol{y}\leftarrow\boldsymbol{x}) < \varepsilon$来决定是否接受采样出来的$\boldsymbol{y}$. Putting all this together, we arrive at the following MCMC sampling procedure:
Metropolis Sampling
The initial state is $\boldsymbol{x}_0$, and the state at time $t$ is $\boldsymbol{x}_t$.
$\boldsymbol{x}_{t+1}$ is sampled via the following procedure:
1. Sample $\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$;
2. Sample $\varepsilon\sim U[0,1]$;
3. Compute $\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)=q(\boldsymbol{x}_t\leftarrow\boldsymbol{y})p(\boldsymbol{y})$;
4. If $\varepsilon \leq \alpha(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$, then $\boldsymbol{x}_{t+1} = \boldsymbol{y}$, otherwise $\boldsymbol{x}_{t+1}=\boldsymbol{x}_t$.
This is the sampling algorithm proposed by Metropolis in 1953, generally referred to as the Metropolis algorithm, or simply the MCMC method. It requires knowing the exact expression for $p(\boldsymbol{x})$, as well as an easy-to-sample transition probability $q(\boldsymbol{y}\leftarrow\boldsymbol{x})$. Then, after iterating the above procedure for enough steps, the resulting $\boldsymbol{x}_t$ can be considered as having been sampled from $p(\boldsymbol{x})$. How many steps count as "enough" here depends on the specific choices of $p(\boldsymbol{x})$ and $q(\boldsymbol{y}\leftarrow\boldsymbol{x})$; generally this can only be determined experimentally, and there's no universal standard.
Some readers might think: isn't this just ordinary rejection sampling? Where does the "piling all the excess probability onto the 'no change' operation" idea show up in practice? In fact, it does differ somewhat from ordinary rejection sampling. In ordinary rejection sampling, if $\varepsilon > \alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})$, we repeat steps 1, 2, 3 over and over until we sample a $\boldsymbol{y}$ that satisfies $\varepsilon \leq \alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})$, and then set $\boldsymbol{x}_{t+1} = \boldsymbol{y}$. In the sampling process above, however, if $\varepsilon > \alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})$, we directly set $\boldsymbol{x}_{t+1} = \boldsymbol{x}_t$, and the value sampled at the next step is $\boldsymbol{x}_{t+2}$. If we assume the distribution becomes stationary starting from time $T$, then $\{\boldsymbol{x}_t\}_{t=T}^{\infty}$ are all samples drawn from $p(\boldsymbol{x})$ (identically distributed but not independent). Clearly, these two different rejection strategies will affect the ultimate distribution of $\{\boldsymbol{x}_t\}_{t=T}^{\infty}$.
It's also worth pointing out that, although the $\boldsymbol{x}$ we deal with in this series of articles are discrete objects, the same conclusions apply equally to sampling continuous objects — we only need to replace the summations in the derivation above with integrals over the probability density; there is no essential difference.
MH Sampling
The Metropolis sampling described above is already usable in many scenarios, but it's not entirely perfect. First, as noted, Metropolis sampling requires knowing the exact expression for $p(\boldsymbol{x})$, which is difficult to obtain in many tasks. Second, the acceptance rate $\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x})$ tends to be too small, leading to long periods of "spinning in place" and thus an excessively long time to reach stationarity. Fortunately, there's a simple trick that solves both problems at once.
The trick is to use the following expression as the acceptance rate:
\begin{equation}\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}) = \min\left(1, \frac{q(\boldsymbol{x}\leftarrow\boldsymbol{y})p(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x})p(\boldsymbol{x})}\right)\end{equation}
This acceptance rate also derives from the identity $\eqref{eq:feihua}$. As far as the acceptance rate is concerned, the only requirement is that it lie between 0 and 1, and beyond that, the closer to 1 the better. To achieve this, we can divide both sides of identity $\eqref{eq:feihua}$ by $\max(\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x}),\alpha(\boldsymbol{x}\leftarrow\boldsymbol{y}))$; that way, one of the two sides becomes 1, while the other, still no greater than 1, gets amplified. Simplifying then gives us $\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x})$. What's rather clever is that $\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x})$ now depends only on the relative value of $p(\boldsymbol{y})$, so we don't need to compute its normalization factor at all.
This gives us an improved version of Metropolis sampling, which is the well-known Metropolis-Hastings Sampling (MH sampling):
Metropolis-Hastings Sampling
The initial state is $\boldsymbol{x}_0$, and the state at time $t$ is $\boldsymbol{x}_t$.
$\boldsymbol{x}_{t+1}$ is sampled via the following procedure:
1. Sample $\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$;
2. Sample $\varepsilon\sim U[0,1]$;
3. Compute $\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t) = \min\left(1, \frac{q(\boldsymbol{x}_t\leftarrow\boldsymbol{y})p(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)p(\boldsymbol{x}_t)}\right)$;
4. If $\varepsilon \leq \mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$, then $\boldsymbol{x}_{t+1} = \boldsymbol{y}$, otherwise $\boldsymbol{x}_{t+1}=\boldsymbol{x}_t$.
Some Analysis
In the previous article, we pointed out that sampling directly from $p(\boldsymbol{x})$ is very difficult — even with rejection sampling, if the dimensionality of $\boldsymbol{x}$ is too high, sampling from the approximate distribution $q(\boldsymbol{x})$ tends to yield a very low acceptance rate, making rejection sampling extremely inefficient and thus unusable. So a natural question arises:
In MCMC methods, we also need to sample from a transition probability $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$, and there is likewise an acceptance probability $\mathcal{A}(\boldsymbol{y}\leftarrow \boldsymbol{x})$. So why is MCMC practical, while ordinary rejection sampling is not?
The reason is that, for direct rejection sampling, what we're sampling directly from the approximate distribution $q(\boldsymbol{x})$ is the entire high-dimensional sequence. If $q(\boldsymbol{x})$ differs substantially from the exact $p(\boldsymbol{x})$, the acceptance probability tends to decay exponentially, making it extremely low. For the MCMC method, however, we don't impose many restrictions on the form of $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$, so we can design $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$ appropriately such that the probability distribution is concentrated only on those $\boldsymbol{y}$ that are similar to $\boldsymbol{x}$. In other words, $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$ is nonzero only when $\boldsymbol{y}$ is quite similar to $\boldsymbol{x}$, and is zero otherwise. As a result, samples drawn from $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$ differ only slightly from the input $\boldsymbol{x}$; since the change is small, the acceptance rate tends to be higher, which is what makes this kind of rejection sampling feasible.
So, put simply, the MCMC method transforms "directly generating $\boldsymbol{x}$" into a gradual, iterative process: "starting from $\boldsymbol{x}_0$, repeatedly tweaking and refining it until a $\boldsymbol{x}$ meeting the conditions is generated." This is what makes it effective, and it's consistent with the "one step at a time" idea mentioned in the previous article.
Next, we'll introduce two examples of MH sampling: Gibbs sampling and simulated annealing, both of which embody this idea of incremental, gradual refinement inherent in MCMC.
Gibbs Sampling
Suppose $\boldsymbol{x}=(x_1,x_2,\dots,x_l)$ is a sequence of length $l$. Gibbs sampling adjusts only one element of it at a time. Specifically, Gibbs sampling defines the reference transition probability as (strictly speaking, one should also divide by $l$, but since this constant doesn't change the result, we omit it):
\begin{equation}q(\boldsymbol{x}_{[x_i=y]}\leftarrow \boldsymbol{x}) = p(y|\boldsymbol{x}_{-i})\triangleq\frac{p(x_1,\dots,x_{i-1},y,x_{i+1},\cdots,x_l)}{\sum\limits_y p(x_1,\dots,x_{i-1},y,x_{i+1},\cdots,x_l)}\end{equation}
where $\boldsymbol{x}_{[x_i=y]}$ is the sequence obtained by replacing the $i$-th position of $\boldsymbol{x}$ with $y$. In other words, it constructs $p(y|\boldsymbol{x}_{-i})$ as the transition probability based directly on the target distribution $p(\boldsymbol{x})$ itself: each time, it first randomly selects a position $i$ from $1,2,\cdots,l$, then replaces the element at position $i$ with a value sampled from the conditional distribution $p(y|\boldsymbol{x}_{-i})$. Even better, we can show that in this case the acceptance probability is always 1:
\begin{equation}\begin{aligned} \mathcal{A}(\boldsymbol{x}_{[x_i=y]}\leftarrow \boldsymbol{x}) =&\, \min\left(1, \frac{p(x_i|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{[x_i=y]})}{p(y|\boldsymbol{x}_{-i})p(\boldsymbol{x})}\right) \\ =&\, \min\left(1, \frac{p(x_i|\boldsymbol{x}_{-i})p(y|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{-i})}{p(y|\boldsymbol{x}_{-i})p(x_i|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{-i})}\right)\\ =&\, \min\left(1, 1\right)\\ =&\,1 \end{aligned}\end{equation}
From this, we can construct the following Gibbs sampling procedure:
Gibbs Sampling
The initial state is $\boldsymbol{x}_0=(x_{0,1},x_{0,2},\cdots,x_{0,l})$, and the state at time $t$ is $\boldsymbol{x}_t=(x_{t,1},x_{t,2},\cdots,x_{t,l})$.
$\boldsymbol{x}_{t+1}$ is sampled via the following procedure:
1. Uniformly sample a position $i$ from $1,2,\cdots,l$;
2. Compute $p(y|\boldsymbol{x}_{t,-i})=\frac{p(x_{t,1},\dots,x_{t,i-1},y,x_{t,i+1},\cdots,x_{t,l})}{\sum\limits_y p(x_{t,1},\dots,x_{t,i-1},y,x_{t,i+1},\cdots,x_{t,l})}$;
3. Sample $y\sim p(y|\boldsymbol{x}_{t,-i})$;
4. $\boldsymbol{x}_{t+1} = {\boldsymbol{x}_t}_{[x_{t,i}=y]}$ (i.e., take $\boldsymbol{x}_{t+1}$ to be the sequence obtained by replacing the $i$-th position of $\boldsymbol{x}_t$ with $y$).
Simulated Annealing
Another example is the Simulated Annealing algorithm. In the previous article, we mentioned that maximizing an objective can actually be viewed as randomly sampling from that objective, and the result of doing so corresponds to the simulated annealing algorithm.
First, suppose the function we want to maximize is $f(\boldsymbol{x})$, and suppose there exists a constant $T > 0$ such that for any $0 < \tau \leq T$, we have
\begin{equation}\sum_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau} < \infty\end{equation}
Here $\tau$ has the physical meaning of temperature. This assumption may look like an extra restriction, but in fact simulated annealing is only applicable to scenarios that satisfy this condition. Since this condition holds, we can construct the distribution
\begin{equation}p_{\tau}(\boldsymbol{x}) = \frac{e^{f(\boldsymbol{x})/\tau}}{\sum\limits_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau}}\end{equation}
Assuming the maximum point is unique, then as $\tau\to 0$, $p_{\tau}(\boldsymbol{x})$ degenerates into a one-hot distribution — that is, only the maximum point has probability 1. If we then sample from $p_{\tau}(\boldsymbol{x})$, the result will be the maximum point. Even when $\tau > 0$, the maximum point still has the highest probability, so sampling from $p_{\tau}(\boldsymbol{x})$ will still tend toward the maximum point. So we can first fix the temperature $\tau$, construct a stochastic process performing MH sampling from $p_{\tau}(\boldsymbol{x})$, and then gradually lower $\tau$ — the final convergence result will then be near the maximum point. This is what's known as "simulated annealing."
To perform MH sampling, we need to construct a transition matrix. The choice used in simulated annealing is fairly simple: based on the current state $\boldsymbol{x}$, design a fixed number of candidates $\boldsymbol{y}$ for the new value of $\boldsymbol{x}$ (this step is usually called "mutation," and it typically just means making a simple modification to $\boldsymbol{x}$ to obtain $\boldsymbol{y}$); then $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$ is chosen uniformly at random from among the candidates. Since it's uniform, $q(\boldsymbol{y}\leftarrow \boldsymbol{x})$ is a constant, and $q(\boldsymbol{y}\leftarrow \boldsymbol{x})=q(\boldsymbol{x}\leftarrow \boldsymbol{y})$, so we have
\begin{equation}\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}) = \min\left(1, \frac{q(\boldsymbol{x}\leftarrow\boldsymbol{y})p_{\tau}(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x})p_{\tau}(\boldsymbol{x})}\right) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x})]/\tau}\right)\end{equation}
So, simulated annealing is a search strategy: if $f(\boldsymbol{y}) \geq f(\boldsymbol{x})$, we update $\boldsymbol{x}=\boldsymbol{y}$; even if not, there's still some probability of updating $\boldsymbol{x}=\boldsymbol{y}$. Throughout the search process, we slowly anneal ($\tau$ gradually approaches 0), and under a suitable annealing schedule, simulated annealing will almost always find the maximum point. Of course, how to design an "appropriate" annealing schedule for a specific problem is itself a question worth pondering.
Simulated Annealing
The initial state is $\boldsymbol{x}_0$, the initial temperature is $\tau_0$, and the temperature decreases according to a predetermined schedule. The state at time $t$ is $\boldsymbol{x}_t$, with temperature $\tau_t$.
$\boldsymbol{x}_{t+1}$ is sampled via the following procedure:
1. Sample $\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$;
2. Sample $\varepsilon\sim U[0,1]$;
3. Compute $\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x}_t)]/\tau_t}\right)$;
4. If $\varepsilon\leq \mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)$, then $\boldsymbol{x}_{t+1} = \boldsymbol{y}$, otherwise $\boldsymbol{x}_{t+1}=\boldsymbol{x}_t$.
If we make a simple modification, changing the acceptance strategy to "accept if $f(\boldsymbol{y}) \geq f(\boldsymbol{x}_t)$, otherwise reject," we get the more basic "Hill Climbing" method. Clearly, hill climbing has a more direct objective and may converge faster early on, but once it falls into a local maximum, it usually can't escape, so its final convergence result isn't as good as simulated annealing's. Of course, one can also improve hill climbing's performance by repeating the experiment with different initial values. Which algorithm to choose ultimately depends on the specific problem.
Summary
This article once again introduces MCMC, a foundational sampling algorithm needed for "constrained text generation." After repeated revisions, I've finally managed to write an introduction to MCMC methods that I'm reasonably satisfied with — the main distinguishing feature being that it answers a number of detail-level questions that are commonly left unaddressed in typical MCMC tutorials. This concludes our introduction to the foundational algorithms. Starting from the next article, we will gradually introduce how to apply these seemingly dry sampling algorithms to vivid, concrete text generation tasks.
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.