起源：Boltzmann神經網絡
Boltzmann神經網絡的結構是由Hopfield遞歸神經網絡改良過來的，Hopfield中引入了統計物理學的能量函數的概念。
即，cost函數由統計物理學的能量函數給出，隨着網絡的訓練，能量函數會逐漸變小。
可視為一動力系統，其能量函數的極小值對應系統的穩定平衡點。
Hinton發明的Boltzmann中乘熱打鐵，對神經元輸出引入了隨機概率重構的概念。其想法來自於模擬退火算法：
首先在高溫下進行搜索，由於此時各狀態出現概率相差不大，系統可以很快進入“熱平衡狀態”，是一種快速找到系統概率的低能區的“粗搜索”過程。
隨着溫度逐漸降低，各狀態出現概率的差距逐漸被擴大，搜索精度不斷提高。最后以一較高置信度達到網絡能量函數的全局最小點。
隨機重構則是Metropolis抽樣提出的，作為模擬退火算法的核心：
也就是說，每次重構，都往Δ能量差減少的方向移動，這是RBM的cost函數設計核心，AutoEncoder也是受此啟發。
Part I：RBM的結構與目標
RBM簡化了Boltzmann網絡，取消了每層神經元之間的連接，即不再是一個遞歸神經網絡，更像是一個MLP。
定義聯合能量函數： $E(v,h)=-b'v-c'h-h'Wv$ (b'、c'分別顯隱bias做矩陣乘法,h'則是h的全部取值，RBM限定了取值范圍{0,1}，由激活函數依據概率分布生成)
由此得到v的邊緣概率分布：$P(v)=\sum_{h}\frac{e^{-E(v,h)}}{Z} \quad\quad Z=\sum_{v}\sum_{h}e^{-E(v,h)}$
Z是歸一化常數，說白了就是保證概率在[0,1]之間湊出的式子，計算cost函數時可以無視。
RBM的目的是算出盡量符合v（即輸入)的概率分布$P(v)$,即實現 $arg\max \limits_{W}\prod_{v\epsilon V } P(v)$
但是問題在於沒有標簽，沒有誤差，無法訓練W，所以無法訓練出$P(v)$的概率分布。
所以早期的RBM采用從h重構v'來計算誤差。重構v'，說的好像挺簡單，但是需要知道$P(v,h)$的聯合概率分布，用這概率分布去生成v'。
通常情況下，聯合概率分布無法計算，但是可以通過條件概率分布+迭代法近似模擬出符合聯合概率的樣本。
RBM中采用的是Gibbs采樣來重構，是MCMC（Monte Carlo Markov Chain）的一種算法特例。
Monte Carlo是根據概率隨機生成數據的算法統稱，Markov Chain則是符合馬爾可夫鏈的算法統稱。【參考】
這樣，只要通過又臭又長的馬爾可夫鏈，使用Gibbs采樣重構出新樣本v'就好了。
有cost函數：$cost=arg\max \limits_{W}\sum_{v\epsilon V}logP(v)-arg\max \limits_{W}\sum_{v'\epsilon V'}logP(v')$
這個cost函數在實際計算時候簡直是折磨，原因有二: ①Z毫無意義且計算麻煩 ② logP(v)牽扯到E(v,h)，E(v,h)的計算也麻煩。
所以Bengio在Learning deep architectures for AI一文中提出了簡化的cost函數——FreeEnergy(自由能量)函數。
定義自由能量函數（無視Z+負對數似然）：
$ \begin{matrix}FreeEnergy(v)=-log(P(v)\cdot Z)\\\\\qquad\qquad\qquad\qquad=-log\sum_{h}e^{-E(v,h)}\end{matrix} $
重定義E函數：$Energy(v)=-b'v-\sum_{i}\sum_{h_{i}}h_{i}(W\cdot v+c),\quad h_{i}\epsilon {\{0,1}\}$
這里的精簡比較神奇，$c'h$被刻意X掉了，論文里假設是h只有一個神經元，這樣可有可無。c的梯度可以在后半部分求出，所以無須擔心c無法更新。
盡管此時的cost函數已經崩壞（Z因子，c'項都不見了），但是骨子還在，對於Gradient Descent影響不大。
這樣，參數對稱的兩部分，可以得到簡化P(v) 【推導詳見論文】:$P(v)=\frac{e^{b'v}}{Z}\prod_{i}\sum_{h_{i}}e^{h_{i}(W\cdot v+c)}$
於是有更加簡化的自由能量函數：$FreeEnergy(v)=-b'v-\sum_{i}log\sum_{h_{i}}e^{h_{i}(W\cdot v+c)},\quad h_{i}\epsilon {\{0,1}\}$
由於取值非0即1，再次簡化, $FreeEnergy(v)=-b'v-\sum_{i}log(1+e^{(W \cdot v+c)})$
最終，使用FreeEnergy(v)替代logP(v)，有簡化cost函數：$cost=FreeEnergy(v))-FreeEnergy(v')$
Part II：重構與Gibbs采樣
由於V層和H層，每層內各個神經元之間並無連接，因而這些神經元的概率分布是獨立的，那么就可以加權通過激活函數計算單個結點的激活概率。
$P(h_{j}=1|v)=Sigmoid(W_{j}\cdot v+c_{j})\\P(v_{i}=1|h)=Sigmoid(W_{i}^{T}\cdot h+b_{i})$
注意反向重構時，同AutoEncoder一樣，仍然使用Tied Weight。
這樣，就有了條件概率分布，可以構建馬爾可夫鏈了。
上圖是一條波動的鏈，v0->h0普通的正向傳播，忽略不計。
從h0，正式開始Gibbs采樣，一個step過程為，hn->vn+1->hn+1，即hvh過程。
$v^{n+1}=Bernoulli-Sigmoid(W^{T}\cdot h^{n}+b)\\h^{n+1}=Bernoulli-Sigmoid(W\cdot v^{n+1}+c)$
當$t \to \infty$時，有$v^{t}\approx v'$
Part III K步對比散度算法（CD-K）、可持久化對比散度算法（Persistent CD)
仔細想想RBM中的循環
①外循環：主迭代，訓練W，近似逼近P(v|h)、P(h|v)。
②內循環，Gibbs采樣，近似逼近P（v）
二重壓力下，RBM網絡的計算量實在可怕。
Hinton在2002年的 Training Products of Experts by Minimizing Contrastive Divergence 一文中有這個式子：
$-\frac{\partial}{\partial W_{ij}}(Q_{0}||Q_{\infty}-Q_{1}||Q_{\infty})\approx<s_{i},s_{j}>Q_{0}-<s_{i},s_{j}>Q_{1}$
即在求v、v'對比相差的cost函數的梯度，用以更新參數時，無限步Gibbs采樣的效果與單步Gibbs采樣效果差距不大。
數學上很難給出解釋，疑似是Gradient Descent與Gibbs Sampling（Gibbs可以看作是GD的先行？）的迭代效果重合。通常取K=1或K=10（更精確）來減少計算量。
既然Gibbs Sampling可以看作是Gradient Descent的先行使者，那么每次做梯度法時候，不必重啟馬爾可夫鏈。
而是利用舊鏈尾產生的hk，作為新鏈頭h0。
數學上仍然很難給出解釋，直覺上說，Gradient Descent的小量更新並不會影響馬爾可夫鏈繼續擴展，即便狗尾續貂，效果仍然不錯。
通過測試，可持久化CD在每次鏈長較長時效果比較好（k=15），若k=1，則影響較大，不建議使用，選擇大量梯度比較穩妥。
Part IV 代理似然函數（Proxies to Likelihood）
Bengio的FreeEnerey方法一定程度上減輕了梯度法的計算量，但是為迭代收斂的評估帶來巨大麻煩。
由於歸一化因子Z的省略，若使用FreeEnerey的cost函數做評估，那么不同數據集的收斂值幾乎是千奇百怪的。
所以需要建立一套數據集公用評估標准，即引入代理似然函數。
代理一：交錯熵（or 最小二乘法）
AutoEncoder中使用的cost函數的修改版，Z變成了馬爾可夫鏈尾的h。
$cost=v\cdot log(Sigmoid(W\cdot h^{-1}+b))+(1-v)\cdot log(1-Sigmoid(W\cdot h^{-1}+b))$
代理二：偽造似然函數（Pseudo likelihood)
仔細回想一下可持久化CD，馬爾可夫鏈狗尾續貂，徹夜未眠。
這樣，若一直使用鏈尾的h做交錯熵，那么這個交錯熵勢必是不停波動的，並不很直觀的。
因而，引入了偽造似然函數：每次選擇v的某一維度值，取反，偽造成重構的v'。
有對數似然函數：$PL(x)=\sum_{i}logP(x_{i}|x_{-i})$
當然並沒有必要那么精確，每次只需隨機抽一個維度進行偽造，乘上N倍即可：
$PL(x)=N\cdot logP(x_{i}|x_{-i}) \quad\quad i\sim Random(0,N)$
$logP(x_{i}|x_{-i})$的條件概率計算頗為麻煩，利用條件概率公式+貝葉斯假設:
$ logP(x_{i}|x_{-i})\approx log\frac{P(x_{i})}{P(x_{i})+P(x_{-i})}\\\\=log\frac{e^{-FreeEnergy(x_{i})}}{e^{-FreeEnergy(x_{i})}+e^{-FreeEnergy(x_{-i})}}\\\\\approx log (Sigmoid(FreeEnergy(x_{-i})-FreeEnergy(x_{i}))) $
Part V 與AutoEncoder的關系
准確來說，AutoEncoder是RBM的簡化衍生物。RBM是一個概率生成模型，而AutoEncoder只是一個普通的模型。
神經網絡的本質是訓練岀能夠模擬輸入的W，這樣，在測試的時候，遇到近似的輸入，W能夠做出漂亮的響應。
RBM選擇概率，是因為有概率論的公式支持。這樣優化網絡，能夠達到上述目標。
只是原始目標不好優化，Hinton才提出對比訓練的方法，即繞了個彎子的選擇重構。
能量函數使得W朝更大概率方向優化。但是，正如線性回歸有最小二乘法和高斯分布兩種解釋一樣。
其實，W的訓練大可不必拘泥於概率，AutoEncoder則繞過了這點，直接選擇了加權重構，所以cost函數簡單。
可以這么說，重構的數學理論基礎就是RBM的原始目標函數。而概率重構啟發了直接重構。兩者近似等價。
從馬爾可夫鏈上看，AutoEncoder可看作是鏈長為1的特殊形式，即一次重構，而RBM是多次重構。
能使用直接重構的另一個原因是，Hinton在實驗中發現，梯度法的一次重構效果出奇好。
所以AutoEncoder中摒棄了麻煩的Gibbs采樣過程。
從GPU計算來看，k=1情況下，AutoEncoer的GPU利用率高(70%),RBM利用率低（30%），一開始實現的時候嚇了一跳。
CUDA執行馬爾可夫鏈效率並不高，目測二項分布隨機重構是由CPU執行的（ _(:_」∠)_想想好像GPU沒那么高級 ）
尤其在把batch_size設為600之后，RBM的GPU利用率居然只有（10%), 所以官方教程把batch_size設為了20，來減小概率生成的計算壓力。
當然k=15時，GPU加速之后仍然十分緩慢。RBM不愧是硬件殺手。
（本圖來自MSI Afterburner，GTX 765M，OC（847/2512/913））
Part VI Gaussian-Bernoulli RBM
仔細看看Bernoulli-Bernoulli的重構過程：
[0/1]->[0/1]->[0/1].....，由於二項分布是離散概率分布，所以重構生成值要么是0，要么是1，似乎有些不靠譜。
所以，就有了Gaussian-Bernoulli RBM，在v層中強制引入連續數值型Gaussian噪聲，主要對原RBM做以下修改：
$E(v,h)=\sum_{i}\frac{(v_{i}-b{i})^{2}}{2\sigma{i}^{2}}-c'h-h'W\frac{v}{\sigma }
\\P(v_{i}=v|h)=N(b_{i}+\sigma _{i}(W\cdot h),\sigma _{i}^{2})
\\P(h_{i}=1|v)=Sigmoid(W\cdot \frac{v}{\sigma }+c)$
其中$Size(\sigma)=NVisble$，為正態分布標准差，為正值，一般不作訓練，取定值1，原因有二：
其一：梯度法訓練會產生負值和0，導致訓練失敗
其二：Gaussian過程並不需要精確的范圍，取的只是正態分布的連續特征，提供Gaussian噪聲。
當然，還是有人提出了可訓練$\sigma$的改進方法，個人實際測試並沒有成功，$\sigma$在梯度下降過程中，出現了過大數值問題，導致網絡崩潰。參考
Gaussian過程極易產生誤解，由於正態分布的取值范圍是R，所以對比散度樣本不可以再使用由Gaussian分布生成v’
Gaussian過程主要為從二項[0/1]的h，提供一個連續型重構緩沖范圍，而不僅僅把馬爾可夫鏈的運行局限在0/1之間
這樣做的原因是，一般輸入是連續型值，符合的是0~1之間的連續概率分布，而不是簡單的離散二項分布。
然后再把Gaussian分布映射回Bernoulli分布h'，最后使用的是h[-1]，令：$v'=Bernoulli-Sigmoid(W\cdot h[-1]+b)$
對比圖如下：
左圖（Bernoulli-Bernoulli） 右圖（Gaussian-Bernoulli)
可以看到，Gaussian分布使得參數學習到了更多的特征，但是也帶來了更大的重構誤差，以及更快下降到局部最小值。學習率和Bernoulli-Bernoulli相同情況下，epoch2就過擬合了。
Hinton在 A practical guide to training restricted Boltzmann machines 提到，Gaussian-Bernoulli模型不穩定，由於Gaussian分布的上界不定，
所以，學習率較之於Bernoulli-Bernoulli應該下降1~2個數量級。同時，應當慎用CD-10、Persistent CD，這些會帶來更不穩定的重構。（個人測試情況，似然幾乎無法下降）