Metatron Dev. X: MCMC

Metropolis-Hastings
#

令 $\mathbf{x}_k$ 为顶点数为 $k$ 的路径, $p_i$ 为 $X_i$ 的密度, 转移函数 $K(\mathbf{x} \to \mathbf{y})$ 为 $\mathbf{x}$ 转移至 $\mathbf{y}$ 的概率密度, 满足 $\sum_{k=1}^\infty\int_{\mathbf{y}_k} K(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k = 1$ , Markov链的演化如下:

\begin{equation} p_i(\mathbf{x}) = \sum_{k=1}^\infty\int_{\mathbf{y}_k} K(\mathbf{y}_k \to \mathbf{x})p_{i-1}(\mathbf{y}_k)\mathrm{d}\mathbf{y}_k \end{equation}

稳态分布(stationary distribution) $p_\infty$ 为不动点即 $p_i = p_{i-1}$ . Metropolis-Hastings以提议密度 $T(\mathbf{x} \to \mathbf{y})$ 生成候选, 以接受概率 $a(\mathbf{x} \to \mathbf{y})$ 决定是否接受:

\begin{equation} p_i(\mathbf{x}) = p_{i-1}(\mathbf{x})(1 - \sum_{k=1}^\infty\int_{\mathbf{y}_k} T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k) + \sum_{k=1}^\infty\int_{\mathbf{y}_k} p_{i-1}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x})\mathrm{d}\mathbf{y}_k \end{equation}

给定目标分布 $\hat{p}$ , 细致平衡(detailed balance)条件如下:

\begin{equation} \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})a(\mathbf{x} \to \mathbf{y}) = \hat{p}(\mathbf{y})T(\mathbf{y} \to \mathbf{x})a(\mathbf{y} \to \mathbf{x}) \end{equation}

归一化常数 $\|\hat{p}\| = \sum_{k=1}^\infty\int_{\mathbf{x}_k} \hat{p}(\mathbf{x}_k)\,\mathrm{d}\mathbf{x}_k$ , 设 $p_{i-1} = \frac{\hat{p}}{\|\hat{p}\|}$ , 代入演化得 $\frac{\hat{p}}{\|\hat{p}\|}$ 为不动点:

\begin{equation} \begin{aligned} p_i(\mathbf{x}) &= p_{i-1}(\mathbf{x}) - \frac{1}{\|\hat{p}\|}\sum_{k=1}^\infty\int_{\mathbf{y}_k}( \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k) - \hat{p}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x}) )\mathrm{d}\mathbf{y}_k\\ &= p_{i-1}(\mathbf{x}) \end{aligned} \end{equation}

按如下方式定义 $a$ , 使得较大的 $\hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})$ 总是被接受:

\begin{equation} a(\mathbf{x} \to \mathbf{y}) = \min(1, \frac{\hat{p}(\mathbf{y})T(\mathbf{y} \to \mathbf{x})}{\hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})}) \end{equation}

不动点不蕴含唯一性, 例如取 $\Omega = \{1, 2, 3\}$ , 转移矩阵行为起点列为终点:

K = \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & 0\\ \frac{1}{2} & \frac{1}{2} & 0\\ 0 & 0 & 1 \end{bmatrix}

令 $p_0 = (\alpha, \beta, \gamma)$ , 一步后即为不动点, 可见不动点与初始状态相关:

p_1 = \left(\frac{\alpha + \beta}{2}, \frac{\alpha + \beta}{2}, \gamma\right)

从任意 $\mathbf{x}$ 出发均能在有限步内到达任何 $\hat{p} > 0$ 的区域则不动点唯一. 若满足细致平衡, 已知 $p_0 = \frac{\hat{p}}{\|\hat{p}\|}$ 不动点为 $\frac{\hat{p}}{\|\hat{p}\|}$ , 因此任意初值均收敛到它:

\begin{equation} \lim_{i \to \infty}p_i = \frac{\hat{p}}{\|\hat{p}\|} \end{equation}

Metropolis Light Transport
#

基于BDPT实现MLT, 令 $p_0$ 为BDPT的PDF, 初始权重 $W_0 = \frac{\hat{p}(X_0)}{p_0(X_0)}$ , 变异保持 $W_i = W_{i-1}$ , 记 $\rho_i(w, \mathbf{x})$ 为 $(W_i, X_i)$ 的联合密度, 加权均衡条件(weighted equilibrium condition)如下:

\begin{equation} \int_\mathbb{R} w\,\rho_i(w, \mathbf{x})\,\mathrm{d}w = \hat{p}(\mathbf{x}) \end{equation}

$i = 0$ 时 $W_0$ 由 $X_0$ 确定, $\rho_0(w, \mathbf{x}) = \delta(w - \frac{\hat{p}(\mathbf{x})}{p_0(\mathbf{x})})p_0(\mathbf{x})$ , 代入直接满足. 变异与 $W$ 无关, 故 $\rho_i$ 与 $p_i$ 的演化相同:

\begin{equation} \rho_i(w, \mathbf{x}) = \rho_{i-1}(w, \mathbf{x})(1 - \sum_{k=1}^\infty\int_{\mathbf{y}_k} T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k) + \sum_{k=1}^\infty\int_{\mathbf{y}_k} \rho_{i-1}(w, \mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x})\mathrm{d}\mathbf{y}_k \end{equation}

两侧乘 $w$ 并对 $w$ 积分, 基于加权均衡条件与细致平衡可得:

\begin{equation} \begin{aligned} \int_\mathbb{R} w\,\rho_i(w, \mathbf{x})\,\mathrm{d}w &= \hat{p}(\mathbf{x}) - \sum_{k=1}^\infty\int_{\mathbf{y}_k}( \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k) - \hat{p}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x}) )\mathrm{d}\mathbf{y}_k\\ &= \hat{p}(\mathbf{x}) \end{aligned} \end{equation}

令 $h$ 为滤波器权重, $I_j$ 为像素 $j$ 的真值, 可验证无偏:

\begin{equation} \begin{aligned} E[W_i h_j(X_i)] &= \sum_{k=1}^\infty\int_{\mathbf{x}_k}\int_\mathbb{R} w\,h_j(\mathbf{x}_k)\rho_i(w, \mathbf{x}_k)\,\mathrm{d}w\,\mathrm{d}\mathbf{x}_k\\ &= \sum_{k=1}^\infty\int_{\mathbf{x}_k} h_j(\mathbf{x}_k)\hat{p}(\mathbf{x}_k)\,\mathrm{d}\mathbf{x}_k\\ &= I_j \end{aligned} \end{equation}

Metropolis-Hastings #

Metropolis Light Transport #

Metropolis-Hastings
#

Metropolis Light Transport
#