Efficient Inference in Fully Connected CRFs with Gaussian Edge Potentials
Intro
回顾一下经典，文章是早期像素级分类做分割效果不好的情况下，利用条件随机场建模，解决细节分割问题的一篇文章，文章主要贡献是提出了Mean Field Approximation来近似原来的全连接条件随机场，同时提出一个Efficient算法将整个更新过程变为线性时间。
Why use CRF
在早期分割领域，像素级分类直接分割的效果并不好，当时并没有太多的结构解决这一问题，而效果不好的原因基本上来自于下采样和网络深度过深导致的信息损失，输出层难以利用这样的信息拟合这样的任务，文章给了一个图来描述这样的情况，其中图(b)就是直接分割的效果：
The Fully Connected CRF Model
为了解决直接分割效果不好的问题，一般利用CRF对图像进行建模，图(c)是基于patch的CRF，概率图中依赖关系只和一定patch内的pixel有关，图(d)是全连接CRF，即一个节点和图中其他所有节点有关。
定义两随机场\(\mathbf{X} = \{X_1, X_2, \cdots X_N \}、\mathbf{I} = \{ I_1, I_2, \cdots I_N \}\)，其中X表示标签变量，I表示图像信息，\(X_i\)定义在\(\mathcal{L} = \{l_1, l_2, \cdots, l_k\}\)，\(I_i\)则表示索引为i的像素所携带的颜色信息。
对于条件随机场\((\mathbf{X},\mathbf{I})\)，其概率分布定义为：
\[P(\mathbf{X}|\mathbf{I}) = \frac{1}{Z(\mathbf{I})} exp(-E(\mathbf{X}, \mathbf{I}))
\]
其中E为能量函数，Z为归一化项，分别定义如下：
\[E(\mathbf{x}|\mathbf{I})=\sum_{i} \psi_{u}\left(x_{i}|\mathbf{I}\right)+\sum_{i<j} \psi_{p}\left(x_{i}, x_{j} | \mathbf{I}\right)
\]
\[Z(\mathbf{I}) = \sum_{x} exp(-E(\mathbf{X}, \mathbf{I}))
\]
其能量函数第一项为一元势函数，描述的是某一节点的势能，可以用神经网络的输出层进行拟合；第二项为二元势函数，描述的是节点i和无向图中的其他所有和i有边的节点的关系势能，而这一项的定义如下：
\[\psi_{p}\left(x_{i}, x_{j}\right| \mathbf{I})=\mu\left(x_{i}, x_{j}\right) \underbrace{\sum_{m=1}^{K} w^{(m)} k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right)}_{k\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right)}
\]
其中\(k\left(\mathbf{f}*{i}, \mathbf{f}*{j}\right)=w^{(1)} \exp \left(-\frac{\left|p_{i}-p_{j}\right|^{2}}{2 \theta_{\alpha}^{2}}-\frac{\left|I_{i}-I_{j}\right|^{2}}{2 \theta_{\beta}^{2}}\right)+w^{(2)} \exp \left(-\frac{\left|p_{i}-p_{j}\right|^{2}}{2 \theta_{\gamma}^{2}}\right)\)为高斯函数，该公式中p表示位置，第一项描述的是颜色相似性和位置相似性，是基于“距离较近且颜色相似的节点可能属于同一标签类”这一先验，第二项描述的是平滑性；f为人工设计的特征，\(w\)为不同类别的权重，\(\mu\)是一指示函数，在两者标签不同时为1，否则为0，即对于标签相同的情况，并不进行能量传导，而标签不通的情况才进行能量传导. \(\theta\)为可学习的参数(\(w\)其实也是可学习的，但是在实验中发现对\(w^{(2)}\))设置为1效果就挺好了，所以这一参数是超参。
关于为什么一元势函数要用神经网路拟合，而二元势函数要在原图上直接利用“人工设计”的特征来描述，是因为CNN本身感受野有限，网络输出feature的某一pixel位置只和其领域一定范围内的像素有关，无法建模全局关系如全局的二元关系，但描述该位置的信息是足够的；而用“人工设计”的特征来描述二元势函数则是根据颜色先验和平滑先验来的，神经网络难以对这两个先验进行建模描述，而我们却可以很方便的通过图像本身运算得到量化描述（尽管需要优化建模的参数）。
对于总像素点数为N的一张图，其边的数量可以达到\(N^2\)，上式的二元势函数的每个节点则需要进行\(N^2\)的运算，本文提出了mean field来近似fully connected crf。
Mean Field Approximation
对于\(P(\mathbf{X}|\mathbf{I})\)，文章希望求的一个分布\(Q(\mathbf{X}|\mathbf{I})\)来近似\(P(\mathbf{X}|\mathbf{I})\)，文章中是用\(P(\mathbf{X})\)和\(Q(\mathbf{X})\)来表示，省略了后面的条件量，表达的意思是一致的。为了和文章保持一致，后面沿用文章的写法。
其中，对于\(Q(\mathbf{X})\)需要满足\(Q(\mathbf{X})=\prod_{i} Q_{i}\left(X_{i}\right)\)，因此只需要最大化分布P和分布Q之间的KL散度，就可以实现对P的近似了。
对于
\[P(\mathbf{X})=\frac{1}{Z} \tilde{P}(\mathbf{X})=\frac{1}{Z} \exp \left(-\sum_{i} \psi_{u}\left(x_{i}\right)-\sum_{i<j} \psi_{p}\left(x_{i}, x_{j}\right)\right) = \frac{1}{Z} \exp \left(-E(\mathbf{X})\right)
\]
这里原文应该把exp里的正负搞错了，纠正一下。
P和Q的KL散度为：
\[\begin{aligned}\mathbf{D}(Q \| P) &=\sum_{\mathbf{x}} Q(\mathbf{x}) \log \left(\frac{Q(\mathbf{x})}{P(\mathbf{x})}\right) \\&=-\sum_{\mathbf{x}} Q(\mathbf{x}) \log P(\mathbf{x})+\sum_{\mathbf{x}} Q(\mathbf{x}) \log Q(\mathbf{x}) \\&=-\mathbf{E}_{\mathbf{U} \sim Q}[\log P(\mathbf{U})]+\mathbf{E}_{\mathbf{U} \sim Q}[\log Q(\mathbf{U})] \\&=-\mathbf{E}_{\mathbf{U} \sim Q}[\log \tilde{P}(\mathbf{U})]+\mathbf{E}_{\mathbf{U} \sim Q}[\log Z]+\sum_{i} \mathbf{E}_{U_{i} \sim Q}\left[\log Q\left(U_{i}\right)\right] \\&=\mathbf{E}_{\mathbf{U} \sim Q}[E(\mathbf{U})]+\sum_{i} \mathbf{E}_{U_{i} \sim Q_{i}}\left[\log Q_{i}\left(U_{i}\right)\right]+\log Z\end{aligned}
\]
第三行到第四行是进行了代换；第四行到第五行是能量函数代换且Z只和P有关，因此可以从期望中提取出来。
此外，Q的约束条件为：
\[\sum_{x_i} Q_i(x_i) = 1
\]
这是因为边际概率和为1。
根据拉格朗日乘数法，构造拉氏函数
\[L = \mathbf{E}_{\mathbf{U} \sim Q}[E(\mathbf{U})]+\sum_{i} \mathbf{E}_{U_{i} \sim Q_{i}}\left[\log Q_{i}\left(U_{i}\right)\right]+\log Z + \lambda \sum_{x_i} (Q_i(x_i) - 1)
\]
因此：
\[\begin{equation}
\begin{aligned}
\frac{\partial L}{\partial Q_i(x_i)} &= \frac{\partial L}{\partial Q_i(x_i)} \left ( \sum_i \mathbf{E}_{U_i \sim Q_i}\left[ \psi_{u}\left(U_{i}\right) \right ] + \sum_{i<j} \mathbf{E}_{U_i \sim Q_i, U_j \sim Q_j} \left[ \psi_{p}\left(U_{i}, U_{j}\right )\right ]
\right.
\\
&+ \left.\sum_{i} \mathbf{E}_{U_{i} \sim Q_{i}}\left[\log Q_{i}\left(U_{i}\right)\right]+\log Z + \lambda \sum_{x_i} (\log Q_i(x_i) - 1) \right) \\
&= \underline{\psi_u(x_i)} + \underline{\frac{\partial L}{\partial Q_i(x_i)}\left(\frac{1}{2} \sum_i Q_i(U_i)\sum_{i\neq j}\mathbf{E}_{U_j \sim Q_j} \left [ \psi_p(U_i, U_j) \right] \right)} \\
&+ \underline{(\log Q_i(x_i) +1)} + \underline{0} + \underline{\lambda} \\
& =\underline{\psi_u(x_i)} + \underline{\left(\frac{1}{2} \sum_{i\neq j}\mathbf{E}_{U_j \sim Q_j} \left [ \psi_p(x_i, U_j) \right] \right)} + \underline{(\log Q_i(x_i) +1)} + \underline{0} + \underline{\lambda} \\
\end{aligned}
\end{equation}
\]
每一项对应的求导我用下划线标出来了，方便区分；第三行是因为将\(i<j\)换成了\(i\neq j\)，所以加上了一个\(\frac{1}{2}\)，之所以换成不等号，是和前面说到的指示函数对应。
接下来让其导数等于0，求解Q：
\[\psi_u(x_i) + \left(\frac{1}{2} \sum_{i\neq j}\mathbf{E}_{U_j \sim Q_j} \left [ \psi_p(x_i, U_j) \right] \right) + (\log Q_i(x_i) +1) + 0 + \lambda = 0
\]
\[\log Q_{i}\left(x_{i}\right)=-\psi_{u}\left(x_{i}\right)-\frac{1}{2} \sum_{j \neq i} \mathbf{E}_{U_{j} \sim Q_{j}}\left[\psi_{p}\left(x_{i}, U_{j}\right)\right]-1-\lambda
\]
所以：
\[\begin{equation}
\begin{aligned}
Q_{i}\left(x_{i}\right)&=\exp \left\{-\psi_{u}\left(x_{i}\right)-\frac{1}{2} \sum_{j \neq i} \mathbf{E}_{U_{j} \sim Q_{j}}\left[\psi_{p}\left(x_{i}, U_{j}\right)\right]-1-\lambda\right\} \\
& =\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{j \neq i} \mathbf{E}_{U_{j} \sim \mathrm{Q}_{j}}\left[\psi_{p}\left(x_{i}, U_{j}\right)\right]\right\}
\end{aligned}
\end{equation}
\]
之所以把\(\frac{1}{2}\)去掉这里我理解是后面的期望项展开是含有可学习参数\(w\)的，因此前面的系数并不重要；后面的常数项被归到normalize项Z里，是因为Q的求和等于1，最终是需要normalize的，可以直接将其归进去而不影响结果。
将相关的变量带入：
\[\begin{aligned}Q_{i}\left(x_{i}=l\right) &=\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{j \neq i} \mathbf{E}_{U_{j} \sim Q_{j}}\left[\mu\left(l, U_{j}\right) \sum_{m=1}^{K} w^{(m)} k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right)\right]\right\} \\&=\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{m=1}^{K} w^{(m)} \sum_{j \neq i} \mathbf{E}_{U_{j} \sim Q_{j}}\left[\mu\left(l, U_{j}\right) k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right)\right]\right\} \\&=\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{m=1}^{K} w^{(m)} \sum_{j \neq i} \sum_{l^{\prime} \in \mathcal{L}} Q_{j}\left(l^{\prime}\right) \mu\left(l, l^{\prime}\right) k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right)\right\} \\&=\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{l^{\prime} \in \mathcal{L}} \mu\left(l, l^{\prime}\right) \sum_{m=1}^{K} w^{(m)} \sum_{j \neq i} k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right) Q_{j}\left(l^{\prime}\right)\right\}\end{aligned}
\]
很不幸的是，对于\(w^{(1)}、\theta\)的学习，并没有利用梯度来学习，主要是速度上的考虑，可能还不如直接暴力搜索快和效果好；文章是利用grid search来搜索的，即预设一些值，利用暴力法挨个尝试在一个固定的验证集上选择最佳参数。
再贴一下文章中收敛的一个大致效果：
文章中还提到指示函数的梯度学习，前面我们将其认为是不可学习的，将结果根据标签是否一致进行一刀切，显然比较hard，文章后面认为指示函数本身可以作为可学习参数进行学习，而这一块可以利用梯度方法进行学习。
对于训练图像\(\mathcal{I}\)和其对应的分割标注\(\mathcal{T}\)，对于学习参数\(\mu\)，其损失表示为：
\[\begin{aligned}\ell\left(\mu: \mathcal{T}^{(n)}, \mathcal{I}^{(n)}\right) &=\log P\left(\mathbf{X}=\mathcal{T}^{(n)} \mid \mathcal{I}^{(n)}, \mu\right) \\&=-E\left(\mathcal{T}^{(n)} \mid \mathcal{I}^{(n)}, \mu\right)-\log Z\left(\mathcal{I}^{(n)}, \mu\right)\end{aligned}
\]
其导数：
\[\begin{aligned}
\frac{\partial}{\partial \mu_{a, b}} \ell\left(\mu: \mathcal{T}^{(n)}, \mathcal{I}^{(n)}\right) &=\frac{\partial}{\partial \mu_{b, a}} \ell\left(\mu: \mathcal{T}^{(n)}, \mathcal{I}^{(n)}\right) \\
&=-\frac{\partial}{\partial \mu_{a, b}} E\left(\mathcal{T}^{(n)} \mid \mathcal{I}^{(n)}, \mu\right)-\frac{\partial}{\partial \mu_{a, b}} \log Z\left(\mathcal{I}^{(n)}, \mu\right) \\
&=-\sum_{m} \frac{1}{2} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) g(\mu_{a,b})-\frac{1}{Z} \frac{\partial}{\partial \mu_{a, b}} Z\left(\mathcal{I}^{(n)}, \mu\right)
\end{aligned}
\]
其中，
\[g(\mu_{a,b}) = 1_{a=\mathcal{T}_{i}^{(n)}} 1_{b=\mathcal{T}_{j}^{(n)}}+1_{b=\mathcal{T}_{i}^{(n)}} 1_{a=\mathcal{T}_{j}^{(n)}}
\]
\[\begin{aligned}
\frac{1}{Z} \frac{\partial}{\partial \mu_{a, b}} Z\left(\mathcal{I}^{(n)}, \mu\right) &=\frac{1}{Z} \sum_{\mathbf{X}} \frac{\partial}{\partial \mu_{a, b}} \tilde{P}\left(\mathbf{X} \mid \mathcal{I}^{(n)}, \mu\right) \\
&=\frac{1}{Z} \sum_{\mathbf{X}} \frac{\partial}{\partial \mu_{a, b}} \exp \left(-E\left(\mathbf{X} \mid \mathcal{I}^{(n)}, \mu\right)\right) \\
&=-\sum_{\mathbf{X}} \frac{1}{Z} \exp \left(-E\left(\mathbf{X} \mid \mathcal{I}^{(n)}, \mu\right)\right) \frac{\partial}{\partial \mu_{a, b}} E\left(\mathbf{X} \mid \mathcal{I}^{(n)}, \mu\right) \\
&=-\sum_{\mathbf{X}} P(\mathbf{X}) \sum_{m} w^{(m)} \frac{1}{2} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right)\left(1_{a=\mathcal{T}_{i}^{(n)}} 1_{b=\mathcal{T}_{j}^{(n)}}+\right. \left.1_{b=\mathcal{T}_{i}^{(n)}} 1_{a=\mathcal{T}_{j}^{(n)}}\right) \\
&=-\sum_{\mathbf{X}} P(\mathbf{X}) \sum_{m} w^{(m)} \frac{1}{2}\left(\sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) 1_{a=\mathcal{T}_{i}^{(n)}} 1_{b=\mathcal{T}_{j}^{(n)}}+\right. \left.\sum_{j \neq i} k^{(m)}\left(\mathcal{T}_{j}^{(n)}, \mathcal{T}_{i}^{(n)}\right) 1_{a=\mathcal{T}_{j}^{(n)}} 1_{b=\mathcal{T}_{i}^{(n)}}\right)\\
&=-\sum_{\mathbf{X}} P(\mathbf{X}) \sum_{m} w^{(m)} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) 1_{a=\mathcal{T}_{i}^{(n)}} 1_{b=\mathcal{T}_{j}^{(n)}}
\end{aligned}
\]
然后再将P换成我们估计的Q：
\[\begin{aligned}&\frac{1}{Z} \frac{\partial}{\partial \mu_{a, b}} Z\left(\mathcal{I}^{(n)}, \mu\right) \\&\approx \sum_{\mathbf{X}} Q(\mathbf{X}) \sum_{m} w^{(m)} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) 1_{a=X_{i}} 1_{b=X_{j}} \\&=\sum_{m} w^{(m)} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) \sum_{\mathbf{X}} Q\left(\mathbf{X} /\left\{X_{i}, X_{j}\right\}\right) 1_{a=X_{i}} Q_{i}\left(X_{i}\right) 1_{b=X_{j}} Q_{j}\left(X_{j}\right) \\&=\sum_{m} w^{(m)} \sum_{i \neq j} k^{(m)}\left(\mathcal{T}_{i}^{(n)}, \mathcal{T}_{j}^{(n)}\right) Q_{i}(a) Q_{j}(b)\end{aligned}
\]
最后得到
\[\begin{aligned}\frac{\partial}{\partial \mu(a, b)} \ell_{n}\left(\mu: \mathcal{I}^{(n)}, \mathcal{T}^{(n)}\right) \approx \sum_{m} w^{(m)}(&-\sum_{i} \mathcal{T}_{i}^{(n)}(a) \sum_{j \neq i} k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right) \mathcal{T}_{j}^{(n)}(b) \\&\left.+\sum_{i} Q_{i}(a) \sum_{j \neq i} k^{(m)}\left(\mathbf{f}_{i}, \mathbf{f}_{j}\right) Q_{i}(b)\right)\end{aligned}
\]
Algorithm
之前我们得到：
\[Q_{i}\left(x_{i}\right)=\frac{1}{Z_{i}} \exp \left\{-\psi_{u}\left(x_{i}\right)-\sum_{j \neq i} \mathbf{E}_{U_{j} \sim \mathrm{Q}_{j}}\left[\psi_{p}\left(x_{i}, U_{j}\right)\right]\right\}
\]
不考虑\(\mu\)的学习的情况下，公式左右两边都是有Q的，因此可以先初始化一个Q，然后不断迭代代入，直到P和Q的KL散度小于一定阈值。
其算法流程如下：
该算法在Message passing阶段的复杂度是\(O(N^2)\)的，利用mean field approximation之后绛为\(O(N)\).
而复杂度的降低来自于对\(k^{(m)}(f_i, f_j)\)用卷积替代，卷积核的大小并非全局(不然不和原来没区别了)，而是利用截断，设置了一个阈值选择卷积核的大小，即只连接高斯分布附近的节点。
替换之后的算法：
可以看到除了用卷积处理，还加上了下采样和上采样，这里先下采样是为了减少运算，卷积过程等价于低通滤波，低通滤波本身就会过滤掉高频部分，试想如果不进行任何采样，那message passing的过程其实就是过滤高频部分，但是加了下采样会使得计算量降低，这里可能会有个疑问，直接下采样不是有可能造成混叠吗，一般我们应该先低通滤波再去下采样的呀？其实这里用到了Nobel恒等式，也就是先低通滤波再下采样和先下采样再低通滤波是等价的，简单的解释是先低通滤波再去下采样，最终留下来的点其实有一些是没用到的，所以可以事先进行采样选择需要的点，再进行滤波，更为详细的理解可以在The Scientist and Engineer's Guide to Digital Signal Processing找到，大概在p40；而上采样部分仅仅是为了恢复形状。
