申探社：不確定度建模（下）-變分推斷及應用

在（上）篇中，我們討論了什么是不確定度，為什么需要關注不確定度建模，以及不確定度可以怎么用。也從最大似然估計（MLE）到最大后驗概率（MAP），講到了貝葉斯推斷(Bayesian Inference)。而我們希望用來建模不確定度的目標模型是貝葉斯神經網絡（BNN），它是一種用神經網絡來建模似然率，然后進行貝葉斯推斷的方法。（中）篇介紹了如何用蒙特卡洛采樣的方法來進行貝葉斯推斷。而本篇會介紹另外一種方法：變分推斷。此外會講解一個非常容易實現的變分推斷特例：MC Dropout。也會討論在實際應用貝葉斯神經網絡（BNN）中的一些問題。本篇涉及內容主要是概念圖中的這一部分：

一. 變分推斷基本思想

回顧一下，在貝葉斯推斷中，我們主要的一個目標是要計算后驗概率?。困難就在于這個后驗概率的解析解很難求。變分推斷的思路是，用另外一個關于??的分布??去作為的近似。而這個?這個分布是參數化的，參數用??表示。我們通過學習參數??從而讓盡可能的接近，然后就可以用作為的近似解了。

明星模仿者，通過不斷修正自己的發型，臉型，口音等等參數，使得自己和明星的差距越來越小，最終就可以作為明星的“近似解”去上綜藝，做節目了。

這樣，我們的目標就變為找到讓盡可能的接近的??值。也就是最熟悉的機器學習最優化問題。我們只需要定義一個Loss，然后讓??成為神經網絡的參數，然后用各種優化方法訓練這個神經網絡就可以了。Loss很好定，我們要讓盡可能的接近，其中一個選擇就是最小化他們的KL散度就好了。也就是說

根據貝葉斯公式

帶入上式可以得到

根據期望值的計算公式??，我們可以把上面的式子變為

如果我們引入一個專門的名稱ELBO，它等于（注意前面有個負號）

則可以把優化目標寫成

因為??與??無關，也與??無關，是個常數，所以有

所以我們的Loss就等于-ELBO就好了。不過ELBO里是三個期望值，而??如果是個連續分布的話，求期望值還得求積分。這個時候我們可以用蒙特卡洛采樣法來幫我們避免求積分（是的，變分推斷里也有蒙特卡洛法，你中有我我中有你）。此時我們的Loss就可以寫成

其中??采樣自??。

Notes 1：我們可以變換下ELBO的表達式，方便我們看看這個表達式內在的含義。

可以發現后面那一項其實就是對數似然率，越大越好。前面那一項則是q分布和先驗分布的KL散度，希望q分布和先驗分布越近越好。似然率+先驗約束，我們的老朋友。

Notes 2: 回憶一下上面的公式

變換下等式，有

因為??，所以有

回憶一下，在貝葉斯公式里，p(D)這一項叫evidence，而ELBO是p(D)取值的下界。這也是為什么叫ELBO（Evidence Lower Bound)的原因，

二. 用變分推斷訓練貝葉斯神經網絡

有了可計算的Loss，我們的建模方法很明確了：我們需要建一個神經網絡，它的可學習參數是??，中間變量??是根據分布??采樣得到的，與之前一樣神經網絡頂層的輸出是似然率??分布的均值（回歸問題中?通常選高斯分布，分類問題則為伯努利分布），根據這個值可以得到樣本整體的對數似然率??，它也是loss的一部分。而loss的其他部分也是關于??的式子，加起來就是上面的??。計算關系如下圖

其中網絡中從??得到??的部分，是一個隨機采樣，假設我們選擇的是高斯分布的話（當然也可以其他分布）例如：

?，其中??指代了?

那么??就是從這個高斯分布中采樣出來的。而這個采樣的過程，在神經網絡里是不可導的，也就是說梯度反向傳播的時候，沒有辦法根據對?的梯度，得到對的梯度，從而更新的值。我們需要用一個叫重參數化（Reparameterization)的trick。

我們把高斯分布做一個變換：?

這樣我們只需要從??采樣一個數，然后乘以??再加上??。如此一來，在做反向梯度傳導的時候，就可以得到對??和??的梯度，從而更新他們的值了。如果選擇其他的分布，也有對應的重參數化方法。

訓練這個神經網絡，我們就可以得到最好的??使得最接近，也就是用作為的近似解?；仡櫼幌仑惾~斯推斷的四個步驟：

第一步：假設先驗??，對似然率建模?
第二步：計算后驗?
展開后有?
第三步：計算預測值的分布?
第四步：計算y的期望值??；
如果需要，計算y的方差??作為預估值的不確定度

我們完成了第二步，對于第三步和第四步，我們需要再次借助蒙特卡洛采樣的方法（不愧是二十世紀十大算法之一）。從中采樣一系列??，然后帶入到貝葉斯神經網絡中，神經網絡的一系列輸出值的均值，就是最終的預估值，方差則可以代表不確定度。對于這一部分，和（中）篇中介紹的蒙特卡洛采樣是一樣的。

三. 變分推斷的代碼示例

（說明：樣例代碼的原則是能說明算法原理，追求可讀性，所以不會考慮可擴展性，性能等因素。為了兼容用pytorch丹爐的朋友，訓練方式用和pytorch更類似的接口。運行環境為python 3，tf 2.3）

1. 構造一些樣本，用來后面訓練和展現效果。（此處參考了兩篇文章中的樣本產生和部分代碼，鏈接：krasserm.github.io/2019 及zhuanlan.zhihu.com/p/10）

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def f(x, sigma):
return 10 * np.sin(2 * np.pi * (x)) + np.random.randn(*x.shape) * sigma

num_of_samples = 64 # 樣本數
noise = 1.0 # 噪音規模

X = np.linspace(-0.5, 0.5, num_of_samples).reshape(-1, 1)
y_label = f(X, sigma=noise) # 樣本的label
y_truth = f(X, sigma=0.0) # 樣本的真實值

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X, y_truth, label='Ground Truth')
plt.title('Noisy training data and ground truth')
plt.legend();

2. 畫出的圖中，實線是y值的真實分布，“+”號是加上一定噪音后，我們觀測得到的樣本，也是后面訓練用的樣本。

3. 與（中）篇的蒙特卡洛法一樣，先驗分布選擇了一個由兩個高斯分布組成的混合高斯分布

不出意外??選擇的是高斯分布，其中??指代了?

其中?

這里對高斯分布的方差做了一些小設計，據說可以幫助??的收斂

from tensorflow.keras.activations import relu
from tensorflow.keras.optimizers import Adam
import tensorflow as tf

class BNN_VI():

def __init__(self, prior_sigma_1=1.5, prior_sigma_2=0.1, prior_pi=0.5):
# 先驗分布假設的各種參數
self.prior_sigma_1 = prior_sigma_1
self.prior_sigma_2 = prior_sigma_2
self.prior_pi_1 = prior_pi
self.prior_pi_2 = 1.0 - prior_pi

# (w0_mu,w0_rho)是用來采樣w0的高斯分布的參數，其他類似
self.w0_mu, self.b0_mu, self.w0_rho, self.b0_rho = self.init_weights([1, 5])
self.w1_mu, self.b1_mu, self.w1_rho, self.b1_rho = self.init_weights([5, 10])
self.w2_mu, self.b2_mu, self.w2_rho, self.b2_rho = self.init_weights([10, 1])

# 把所有的mu和rho放在一起好管理，模型里可學習參數是mu和rho，不是w和b
self.mu_list = [self.w0_mu, self.b0_mu, self.w1_mu, self.b1_mu, self.w2_mu, self.b2_mu]
self.rho_list = [self.w0_rho, self.b0_rho, self.w1_rho, self.b1_rho, self.w2_rho, self.b2_rho]
self.trainables = self.mu_list + self.rho_list

self.optimizer = Adam(0.08)

def init_weights(self, shape):
# 初始化可學習參數mu和rho們
    w_mu = tf.Variable(tf.random.truncated_normal(shape, mean=0., stddev=1.))
    b_mu = tf.Variable(tf.random.truncated_normal(shape[1:], mean=0., stddev=1.))
    w_rho = tf.Variable(tf.zeros(shape))
    b_rho = tf.Variable(tf.zeros(shape[1:]))
return w_mu, b_mu, w_rho, b_rho

def sample_w_b(self):
# 根據mu和rho們，采樣得到w和b們
self.w0 = self.w0_mu + tf.math.softplus(self.w0_rho) * tf.random.normal(self.w0_mu.shape)
self.b0 = self.b0_mu + tf.math.softplus(self.b0_rho) * tf.random.normal(self.b0_mu.shape)
self.w1 = self.w1_mu + tf.math.softplus(self.w1_rho) * tf.random.normal(self.w1_mu.shape)
self.b1 = self.b1_mu + tf.math.softplus(self.b1_rho) * tf.random.normal(self.b1_mu.shape)
self.w2 = self.w2_mu + tf.math.softplus(self.w2_rho) * tf.random.normal(self.w2_mu.shape)
self.b2 = self.b2_mu + tf.math.softplus(self.b2_rho) * tf.random.normal(self.b2_mu.shape)
self.w_b_list = [self.w0, self.b0, self.w1, self.b1, self.w2, self.b2]

def forward(self, X):
self.sample_w_b()

# 簡單的3層神經網絡結構
    x = relu(tf.matmul(X, self.w0) + self.b0)
    x = relu(tf.matmul(x, self.w1) + self.b1)
self.y_pred = tf.matmul(x, self.w2) + self.b2
return self.y_pred

def prior(self, w):
# 先驗分布假設
return self.prior_pi_1 * self.gaussian_pdf(w, 0.0, self.prior_sigma_1) \
+ self.prior_pi_2 * self.gaussian_pdf(w, 0.0, self.prior_sigma_2)  

def gaussian_pdf(self, x, mu, sigma):
return tf.compat.v1.distributions.Normal(mu,sigma).prob(x)

def get_loss(self, y_label):
self.loss = []
for (w_or_b, mu, rho) in zip(self.w_b_list, self.mu_list, self.rho_list):
# 這里的q_theta_w和文章對應，其中的w指的是(mu,rho), 而w_or_b中的w就是權重中的w
        q_theta_w = tf.math.log(self.gaussian_pdf(w_or_b, mu, tf.math.softplus(rho)) + 1E-30)
        p_theta = tf.math.log(self.prior(w_or_b) + 1E-30)
# 公式中三項中的兩項
self.loss.append(tf.math.reduce_sum(q_theta_w - p_theta))

    p_d_theta = tf.math.reduce_sum(tf.math.log(self.gaussian_pdf(y_label, self.y_pred, 1.0) + 1E-30))
# 公式中三項中的另外一項
self.loss.append(tf.math.reduce_sum(-p_d_theta))
return tf.reduce_sum(self.loss)

def train(self, X, y_label):
    loss_list = []
for _ in range(2000):
with tf.GradientTape() as g:
self.forward(X)
        loss = self.get_loss(y_label)
      gradients = g.gradient(loss, self.trainables)
self.optimizer.apply_gradients(zip(gradients, self.trainables))
      loss_list.append(loss.numpy())
return loss_list

def predict(self, X):
return [self.forward(X) for _ in range(300)]

X = X.astype('float32')
y_label = y_label.astype('float32')

model = BNN_VI()
loss_list = model.train(X,y_label)
plt.plot(np.log(loss_list))
plt.grid()

4. 畫出來的曲線是每一輪的loss

5. 預估階段，我們把預估的均值和不確定度都畫出來

X_test = np.linspace(-1., 1., num_of_samples * 2).reshape(-1, 1).astype('float32')
y_preds = model.predict(X_test)
y_preds = np.concatenate(y_preds, axis=1)

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X_test, np.mean(y_preds, axis=1), 'r-', label='Predict Line')
plt.fill_between(X_test.reshape(-1), np.percentile(y_preds, 2.5, axis=1), np.percentile(y_preds, 97.5, axis=1), color='r', alpha=0.3, label='95% Confidence')
plt.grid()
plt.legend()

6. 圖中，紅色線條就是所有BNN輸出值的均值，也就作為最終的預估值。而紅色的區域寬窄，則反應了所有BNN輸出值的不確定度（為了方便可視化這里用分位數）。可以看到，結果和（中）篇中蒙特卡洛采樣法得到的結果是類似的：對于沒有樣本的區域，不確定度較大，而對于樣本比較密集的地方，不確定度小。另外，在樣本有覆蓋的領域，轉彎的地方，因為要偏離原來的路線，不確定度大；而直線的地方，則不確定度更小。

四. MC Dropout

前面介紹的蒙特卡洛采樣法和變分推斷的方法，實現起來都略顯復雜。對比之下，接下來介紹的這個方法，實現非常簡單，適合工業界應用。

這個方法的做法簡單到只有兩句話：在原有的只預估均值的神經網絡里，為每一層的所有權重都添加Dropout層，然后在預估（inference）的時候，讓dropout繼續生效。這樣同一個樣本，每次inference都會預估出不一樣的值（因為有dropout），把這些值的均值，就作為預估值，這些值的方差，就作為預估值的不確定度。（Paper鏈接：proceedings.mlr.press/v）

這個做法，可以證明效果和變分推斷是等效的。不過盡管做法這么簡單，但是這個證明卻非常復雜，感興趣的朋友可以看這個證明文檔（文檔鏈接：arxiv.org/pdf/1506.0215）

MC Dropout就像在巨無霸漢堡的每一層肉中間，都加上一層生菜。

我們介紹的訓練貝葉斯神經網絡的三種方法的異同如下

對于MC Dropout來說，也可以認為在Inference階段，??是根據為以下分布采樣得到的：

其中??每一緯度互相獨立，且為Bernoulli分布

其中k為dropout的保留率。

因此MC Dropout作為一種變分推斷的一個特例，在這一點上是一致的，?都是從某個以??為參數的分布中采樣出來。

五. MC Dropout代碼示例

1.（這里略去了生產樣本的代碼，和上面普通變分推斷的代碼是精確一致的）我們直接用keras的接口，像普通神經網絡一樣把結構搭建起來，唯一的不同是每層之間都插入了Dropout層。

import tensorflow as tf
from tensorflow.keras.optimizers import Adam

model = tf.keras.models.Sequential([
  tf.keras.layers.Flatten(input_shape=(1,1)),
  tf.keras.layers.Dense(30, activation='tanh'),
  tf.keras.layers.Dropout(0.3),
  tf.keras.layers.Dense(30, activation='tanh'),
  tf.keras.layers.Dropout(0.3),
  tf.keras.layers.Dense(1, activation='linear')
])

model.compile(optimizer=Adam(learning_rate=0.1),
              loss='mean_squared_error',
              metrics=['MeanSquaredError'])

model.fit(X, y_label, epochs=2000)

2. 在預測的時候，通過 training=True 這個設置，使得Dropout層在inference的時候，也會概率丟棄一部分節點。不設置這個參數的話，dropout層在inference的時候，默認是會保留所有的節點，不執行dropout。

X_test = np.linspace(-1., 1., num_of_samples * 2).reshape(-1, 1).astype('float32')
y_preds = [model(X_test, training=True) for _ in range(300)]
y_preds = np.concatenate(y_preds, axis=1)

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X_test, np.mean(y_preds, axis=1), 'r-', label='Predict Line')
plt.fill_between(X_test.reshape(-1), np.percentile(y_preds, 2.5, axis=1), np.percentile(y_preds, 97.5, axis=1), color='r', alpha=0.3, label='95% Confidence')
plt.grid()
plt.legend()

3. 圖中，紅色線條就是所有BNN輸出值的均值，也就作為最終的預估值。而紅色的區域寬窄，則反應了所有BNN輸出值的不確定度（為了方便可視化這里用分位數）。可以看到，結果和（中）篇中蒙特卡洛采樣法或普通變分推斷的結果是一樣的。

六. 不確定度建模實際應用中的一些問題

(1) 對現有系統侵入性

不管是蒙特卡洛采樣法還是普通的變分推斷，對原有的CXR模型來說，都需要對模型結構或者訓練流程做較大的改動。這也是阻攔不確定性建模落地的主要原因。而MC Dropout對原有模型的侵入型很小，比較方便落地。對于添加的dropout層，也許本來也可以幫助模型抑制過擬合，而提高預估的準確性。

(2) 性能問題

相比普通的模型，貝葉斯神經網絡要預估出很多個值，再求方差來作為不確定度的度量，在線inference時的計算量會翻幾倍，對性能會有較大影響?？梢钥紤]做以下幾點的優化來緩解性能問題：

1. 只在最頂上幾層添加dropout做tradeoff。雖然理論上不能保證和變分完全等效，但是在筆者的實驗中，效果也是可接受的。

2. 原有模型不變，仍舊只預估期望值。單獨使用一個結構更簡單的模型，再加上dropout來預估不確定度。可以這么做是因為不確定度的精度不需要那么高。另外這樣的好處是可以完全和原有模型結耦，0侵入性。但是缺點是這個更簡單的模型，和原有模型不一致，不確定度和預估值不是一起訓練出來的，可能會“不配套”。

3. 在線inference時每次inference是獨立的，可以進行并行化。對于變分推斷（包括MC Dropout）可以先把??的采樣結果存儲起來，這樣就可以和蒙特卡洛采樣法一樣，可以把貝葉斯神經網絡看成是多個模型的ensemble。

4. 事實上，在很多實際應用場景，因為對不確定度的預估準確度要求并不高，因此對實時性的要求也不高，我們可以通過離線計算來緩存不確定度，然后定期更新就好。例如我們主要關心新廣告的不確定度，則可以定期采樣一些涉及該廣告的請求，然后離線計算出這些請求的不確定度，再取平均（抹平人群和上下文差異）作為該廣告的不確定度。如果需要關心某人群與新廣告組合的不確定度，也可以如法炮制，離線統計人群與新廣告交叉緯度的平均不確定度。

（3）不確定度的預估準確性評估

貝葉斯神經網絡模型中用了不同的網絡結構，或者不同的超參數，哪個模型預估的不確定度更準確呢？要直接評估不確定度預估的準確性不太容易，我們只好通過兩個間接評估的方法來評估。

直接評估不好評估，則考慮間接評估的例子如下。

方法一：如果我們是通過預估y值的分布，再用方差來作為不確定度的度量，則可以通過評估y值分布的準確性來間接評估方差的準確性。我們先用訓練集進行模型訓練，然后在測試集上做評估。對有n個樣本的測試集中每個樣本??的特征??，根據m個已存儲（當使用MCMC）或者采樣出（當使用變分推斷或MC Dropout）的??的值，經過模型預估出m個?分布，然后計算??在這m個分布的平均log似然率。再對所有樣本求出總體的平均log似然率。即

如果我們的任務是regression，如前文所述，神經網絡的輸出為高斯分布?的均值，方差為常數??。則上式變為

測試集的平均log似然率越高，說明預估的y分布，和實際分布越符合。從直覺來看，當對某個樣本的y值的預估方差較小時，如果實際的y值方差較大，也就是落在離均值比較遠的地方，則會有一個較低的似然率，從而懲罰偏小的預估方差；而當預估方差較大時，如果實際的y值方差較小，即y值都落在均值附近的地方，也會因為均值附近概率密度被分掉一部分到其他地方，而得到一個較低的似然率，從而懲罰偏大的預估方差。

不過平均log似然率評估的是整個分布的預估準確度，也包括了均值的預估準確度。所以log似然率的提升，不一定是不確定度預估準確性的提升。

方法二：直接評估不好評估，我們只好通過使用了不確定度預估值的任務的效果來間接評估了。下面用冷啟動作為一個例子，但是其他的任務也都可以設計相應的方法。在冷啟動中，我們使用不確定度來決定探索哪些廣告，那么我們就可以根據探索的效率來評估不確定度預估的準確度。

我們可以先用第1到第k天的樣本訓練一個Base模型?模型。然后分別執行一個依賴了不確定度預估模型A和模型B的某個探索樣本選擇策略，從第k+1到第m天的樣本中進行探索廣告的挑選，分別根據模型A和模型B選擇出樣本數量相同的樣本集??和??。用這兩個樣本集對?模型進行增量訓練，得到模型??和??。然后用第m+1天到第n天的測試集，來評估這兩個模型??和??以及Base 模型的AUC。如果模型?的AUC比模型的AUC提高的程度大于模型B比Base模型，則說明不確定預估模型A的預估效果可能好于模型B，或者說至少更適合這個探索策略。

上面的這個流程的原理，是不是和Meta Learning的有些類似？如果我們把Base模型的模型結構或超參數，不確定度模型的結構或者超參數，以及探索樣本選擇策略，都作為meta模型的可學習參數，我們就可以試圖去學一個好的學習方法，這個學習方法好不好，則可以通過在最終測試集上的AUC等metrics來評估。最終得到一個較好的學習方法，可以在同樣的探索成本下，獲得最高的模型效果提升，也就是冷啟動探索效率的提升。

當然，實際線上的冷啟動任務比構造的這個任務要更復雜，例如在這個任務中，通過限制探索樣本的數量來限制探索所需要的成本，而實際線上搜集每個樣本的成本是不相同的，數量一致不代表資源一致。

講到這里，本系列的主要內容就介紹完了，感謝大家的支持。希望大家共同來探索不確定度建模在廣告系統中的應用。

還是廣告

又到了招聘旺季，我們正在大力尋找志同道合的朋友一起在某手商業化做最有趣最前沿的廣告算法，初中高級廣告算法職位均有HC（迅速上車，還能趕上上市）。感興趣的朋友歡迎加我的個人微信約飯約咖啡索要JD或發送簡歷。（產品和運營也在招人，看機會的朋友我可以幫忙直接推薦給leader們。）

作者個人微信（添加注明申探社讀者及簡單介紹）：

歡迎掃描下面二維碼關注本公眾號，也歡迎關注知乎“申探社”專欄