前言

网路上有许多介绍VAE的文章与影片，但许多解释公式都无法得知为什么，但这也是个人对于看paper的功力太差所以才能努力爬文观看他人解析，而最后终于找到一个合乎我疑问的解释[1]，最后也发现[2]的解释方法就是利用[1]来去推导，最终解开了我的为甚么。
机率公式1
机率公式2

---

20190706更新:

修正相关叙述。新增先后验机率。

VAE最大化的困难

从[4]得知，假设平常训练模型输入为X，输出一个类别Y，则我们会建立一个模型f(X;W)，使得输入X输出Y的机率最大，即是p(Y|X)最大。
而若想使用生成模型，则必须转为贝叶斯公式，如下图(z=Y)。这里X为观测变量意旨容易观察到X，z为隐藏变量意旨不易观察，于是我们需要知道z的隐藏变量就可使用贝叶斯公式求解。
而p(z)的z机率能由我们决定，而较难解决的则是p(X|z)，这能使用最大似然估计解这个参数。

来源[4]。

但若隐藏变量是一个连续机率分布的空间，公式如下图，分母意思是一个连续机率为隐藏变量z输出X的机率乘上z的机率，而这在数学上要求导是有难度的。

来源[4]。

简单来说求解最大似然估计在一个高维空间里要去解这个公式就不那么容易。而重要的是最后要去度量与p(z|X)的差异，但这里我们并不知道X的真实分布是什么(另一篇文章是用p(x)来解说)。

然而不管是要近似文章一的p(x)或是文章二的p(z|X)只要知道最后我们都必须去度量与真实分布差异，所以迴避计算高维的问题最后VAE当中使用一个机率分布q(z)来逼近p(z|X)或p(x)。

先后验机率

在这里先验机率可以想像成是一个真实存在的，但我们无法得知他真实的机率是多少，而后验机率则是在某个事件当中发生另一个事件的机率是多少，而在VAE当中就能利用这特性去近似函数(实际上要最大似然)。
主要是先验机率是存在的我们无法去控制，若可控制则会改变原先的结果，但后验则可想像成一个类似映射的方式，多添加一个可控制的机率函数产生出原先的真实的机率。而这边使用贝叶斯则可转换为这样的构思。

首先将上述贝叶斯公式转换为文字叙述，使得快速理解。

p(z|X)为在X事件中发生z的机率。因是由X当中取出z，在机率当中可称为z的后验机率。p(X|z)为在z事件中发生X的机率。因是由z当中取出X，在机率当中可称为X的后验机率。p(z)为发生z的机率。在机率当中可称为z的先验机率。p(X)为发生X的机率。在机率当中可称为X的先验机率。

[6]提到后验与先验使用贝叶斯公式则

比例P(B|A)/P(B)也有时被称作标准似然度（standardised likelihood），贝氏定理可表述为：
后验概率 = 标准似然度*先验概率

在上述中可得知虽然很难计算出最大似然数，但可以假设一个q(z)，而p(z|X)是后验机率，则只要使用KL散度最小化q(z)与p(z|x)即可。(只要Encoder输出的z是高斯分布就足够了)

ELBO

首先介绍ELBO，其中有使用到贝叶斯(可看成一般机率公式转换)，这里使用一个q(z)趋近p(z|x)，而使用的度量方式为KL散度，KL散度即是两者分布的差异。但我们无法直接求min，因为我们无法得知p(x)的真实分布情况只知道它是存在的，如下图的推导。

公式，来源[1]。

1.KL散度指数乘法转为减法。
2.p(z|x)转换贝叶斯公式。

推导，来源[1]。

但可将上述公式转换为ELBO进行推导。如下图。

1.p(x)扣掉上述推导后公式，原先的p(x)就消失了，这动作主要为了后面继续推导。

来源[1]。

你会好奇为何可以直接扣除logp(x)，首先观看[1]的说法

由于logp(x)相对于q(z)是一个常量，而我们想要最小化KL，ELBO等于负的KL加上一个常量，所以我们最大化ELBO就等价于最小化KL。最大化logp(x)就是对观测数据的极大似然估计（即log evidence），由于KL是非负的，所以这个目标函数是极大似然估计的下确界（即evidence lower bound, ELBO）。

简单来说我们要最小化q(z)(微分)，而logp(x)它不依赖z(微分是常数)所以在这里可视为常量，那最小化KL即是最大化ELBO，因ELBO=常量-KL，KL越小ELBO越大。接着继续推导化简为如下。

1.p(z,x)机率公式转换为p(x|z)p(z)。
2.指数关係，乘法变加法。
3.后两项转为KL。

来源[1]。

最后这公式没有了计算p(x)分布，而是用一个隐变量z来近似x，而在VAE当中这隐变量都假设为方差为1平均为0。

VAE

ELBO在VAE运用

这里介绍两种方法，方法一来自[2]，方法二则是使用上面公式带入，方法一其实与方法二很像，但方法一有些数学公式没很熟悉则会转不过来。因此将上述介绍的公式直接使用在上面会更加清晰。

方法一[2]

作者使用联合机率分布，以机率来说则在计算同时产生x和z的机率，Z是一个隐藏变量(简单来说就是Encoder后的输出)，q(X,Z)则是一个近似p(X,Z)，最后面则会直接假设Z输出是方差=1、均值=0的二维高斯变量。

假设q(x, z)是一个度量p(x, z)的联合机率分布，使用KL散度进行度量计算，如下图。

来源[2]。

1.将p(x,z)使用机率转换公式转为p(x)p(z|x)。
2.转为期望值公式。
3.将右边公式p(x)乘法转为加法。
4.分开为左右两式。
5.左式直接对Z做微分，可以看到Z为一个输出乘上p(x)所以结果为p(x)。([2]解释为左式因p(x)与dz无关所以提出)。

来源[2]。

[2]解说为，lnp(x)虽然是一个複杂的分布，虽然无法求出但可确定它是存在，所以可当作一常数。

来源[2]。

1.q(x,z)转为q(x|z)q(z)。
2.将右式q(x|z)提出。
3.转为期望公式。
4.右式转为KL。

来源[2]。

方法二

为了验证[1]所以使用上述ELBO求出。假设q(z)近似p(z|X)。

1.展开KL散度。
2.提出q(z)。
3.logp(z|X)使用贝叶斯转换。
4.将与z无关的项提出，logp(x)。

将logp(x)移置左项，最小化KL等于最大化ELBO。

1.带入上述推导公式。
2.logp(z,X)使用机率公式转为logp(X|z)p(z)。
3.提出logp(z)，指数乘法等于加法。
4.将q(z)和p(z)转为KL。

ELBO展开公式

接着要将最后推倒结果给展开，这里延续方法一[2]。

KL项

对于KL项假设q(z)是个正态分布，方差=1、均值=0，这里公式基本上都是KL展开化简，有兴趣可前往详细KL多维高斯分布推导。

来源[2]。

左项

论文使用了两种分布，其第一种结果为交叉熵，第二种为均方误差。

第一种
使用伯努利分布，如下。
1.将每个维度的机率带入伯努利分布。
2.连乘转为指数运算。

来源[2]。

结果即是交叉熵。这里要注意的是ρ(z)是Decoder运算，因伯努利分布所以输出资料必须在0~1，可用sigmoid函数实现。

第二种
使用正态分布，而这里与上述正态分布推导相同。
1.带入正态分布。
2.连乘转为指数运算。
3.固定方差为常数，因固定方差所以后向可无视，而平方中的alpha可提出。

程式码

这里的程式码与AutoEncoder很像，差别在于Encoder和loss不同。

参数

learning_rate: 学习率。
batch_size: 批次训练数量。
train_times: 训练次数。
train_step: 验证步伐。
encoder_hidden: 编码神经网路隐藏层数量。
decoder_hidden: 解码神经网路隐藏层数量。
output_size: 编码输出数量。

learning_rate = 0.001batch_size = 128train_times = 200train_step = 1encoder_hidden1_size = 500encoder_hidden2_size = 300encoder_hidden3_size = 150decoder_hidden1_size = 150decoder_hidden2_size = 300decoder_hidden3_size = 500output_size = 2

神经网路层

批量归一化
之前有介绍，可前往网址观看。

def layer_batch_norm(x, n_out, is_train):    beta = tf.get_variable("beta", [n_out], initializer=tf.zeros_initializer())    gamma = tf.get_variable("gamma", [n_out], initializer=tf.ones_initializer())    batch_mean, batch_var = tf.nn.moments(x, [0], name='moments')    ema = tf.train.ExponentialMovingAverage(decay=0.9)    ema_apply_op = ema.apply([batch_mean, batch_var])    def mean_var_with_update():        with tf.control_dependencies([ema_apply_op]):            return tf.identity(batch_mean), tf.identity(batch_var)    mean, var = tf.cond(is_train, mean_var_with_update, lambda:(batch_mean, batch_var))    x_r = tf.reshape(x, [-1, 1, 1, n_out])    normed = tf.nn.batch_norm_with_global_normalization(x_r, mean, var, beta, gamma, 1e-3, True)    return tf.reshape(normed, [-1, n_out])

神经网路层

def layer(x, weights_shape, activation='relu'):    init = tf.random_normal_initializer(stddev=np.sqrt(2. / weights_shape[0]))    weights = tf.get_variable(name="weights", shape=weights_shape, initializer=init)    biases = tf.get_variable(name="biases", shape=weights_shape[1], initializer=init)    mat_add = tf.matmul(x, weights) + biases    #mat_add = layer_batch_norm(mat_add, weights_shape[1], tf.constant(True, dtype=tf.bool))    if activation == 'relu':        output = tf.nn.relu(mat_add)    elif activation == 'sigmoid':        output = tf.nn.sigmoid(mat_add)    elif activation == 'softplus':        output = tf.nn.softplus(mat_add)    else:        output = mat_add    return output

Encoder

这里输出两个向量，一个为标準差，一个为均值，依照原先构思乘上一个正态分布。而很多人也会将标準差先改为ln(std^2)，但这里使用上述推导公式标準差实作。

def encoder(x):    with tf.variable_scope("encoder"):        with tf.variable_scope("hide1"):            hide1 = layer(x, [784, encoder_hidden1_size], activation='softplus')        with tf.variable_scope("hide2"):            hide2 = layer(hide1, [encoder_hidden1_size, encoder_hidden2_size], activation='softplus')        with tf.variable_scope("hide3"):            hide3 = layer(hide2, [encoder_hidden2_size, encoder_hidden3_size], activation='softplus')        with tf.variable_scope("mean"):            mean = layer(hide3, [encoder_hidden3_size, output_size], activation=None)        with tf.variable_scope("std"):            log_var = layer(hide3, [encoder_hidden3_size, output_size], activation=None)        with tf.variable_scope("output"):            eps = tf.random_normal(shape=tf.shape(log_var), mean=0, stddev=1, dtype=tf.float32)            # log_var = ln(std^2)            #output = mean + tf.exp(log_var) * 0.5 * eps            # log_var = std            output = mean + log_var * eps    return output, mean, log_var

Decoder

def decoder(x):    with tf.variable_scope("decoder", reuse=tf.AUTO_REUSE):        with tf.variable_scope("hide1", reuse=tf.AUTO_REUSE):            hide1 = layer(x, [output_size, decoder_hidden1_size], activation='softplus')        with tf.variable_scope("hide2", reuse=tf.AUTO_REUSE):            hide2 = layer(hide1, [decoder_hidden1_size, decoder_hidden2_size], activation='softplus')        with tf.variable_scope("hide3", reuse=tf.AUTO_REUSE):            hide3 = layer(hide2, [decoder_hidden2_size, decoder_hidden3_size], activation='softplus')        with tf.variable_scope("output", reuse=tf.AUTO_REUSE):            output = layer(hide3, [decoder_hidden3_size, 784], activation='sigmoid')    return output

预测函数

def predict(x):    output, mean, log_var = encoder(x)    return decoder(output), mean, log_var

损失函数

依照上述推导公式带入，若再Encoder使用ln(std^2)则要使用注解公式。

def loss(output, x, mean, log_var):    cross_loss = -tf.reduce_sum(x * tf.log(output + 1e-8) + (1 - x) * tf.log(1 - output + 1e-8), 1)    cross_mean_loss = tf.reduce_mean(cross_loss)    # log_var = ln(std^2)    #kl_loss = 0.5 * tf.reduce_sum(tf.square(mean) + tf.exp(log_var) - log_var - 1, axis=1)    # log_var mean std    kl_loss = 0.5 * tf.reduce_sum(tf.square(mean) + tf.square(log_var) - tf.log(tf.square(log_var)) - 1, axis=1)    kl_mean_loss = tf.reduce_mean(kl_loss)    loss = cross_mean_loss + kl_mean_loss    loss_his = tf.summary.scalar("loss", loss)    return loss

训练函数

def train(loss, index):    return tf.train.AdamOptimizer(learning_rate=learning_rate, beta1=0.9, beta2=0.999, epsilon=1e-08).minimize(loss, global_step=index)    #return tf.train.RMSPropOptimizer(learning_rate, decay=0.9).minimize(loss,    #global_step=index)    #return tf.train.AdagradOptimizer(learning_rate).minimize(loss,    #global_step=index)    #return tf.train.MomentumOptimizer(learning_rate,    #momentum=0.9).minimize(loss, global_step=index)    #return tf.train.GradientDescentOptimizer(learning_rate).minimize(loss,    #global_step=index)

验证函数

def image_summary(label, image_data):    reshap_data = tf.reshape(image_data, [-1, 28, 28, 1])    tf.summary.image(label, reshap_data)def accuracy(output, t):    image_summary("input_image", t)    image_summary("output_image", output)    sqrt_loss = tf.sqrt(tf.reduce_sum(tf.square(tf.subtract(output, t)), 1))    y = tf.reduce_mean(sqrt_loss)    tf.summary.scalar("accuracy error", y)    return y

训练

1.初始化资料。
2.预测函数。
3.损失函数。
4.训练函数。
5.验证函数。
6.tensorboard可视化。
7.开始训练。
8.输出q(X|z)结果。

if __name__ == '__main__':    # init    mnist = input_data.read_data_sets("MNIST/", one_hot=True)    input_x = tf.placeholder(tf.float32, shape=[None, 784], name="input_x")    # predict    predict_op, mean, log_var = predict(input_x)    # loss    loss_op = loss(predict_op, input_x, mean, log_var)    # train    index = tf.Variable(0, name="train_time")    train_op = train(loss_op, index)    # accuracy    accuracy_op = accuracy(predict_op, input_x)    # graph    summary_op = tf.summary.merge_all()    session = tf.Session()    summary_writer = tf.summary.FileWriter("log/", graph=session.graph)    data_x = tf.placeholder(tf.float32, shape=[None, output_size], name="data_x")    decoder_op = decoder(data_x)    init_value = tf.global_variables_initializer()    session.run(init_value)    for time in range(train_times):        avg_loss = 0.        total_batch = int(mnist.train.num_examples / batch_size)        for i in range(total_batch):            minibatch_x, _ = mnist.train.next_batch(batch_size)            session.run(train_op, feed_dict={input_x: minibatch_x})            avg_loss += session.run(loss_op, feed_dict={input_x: minibatch_x}) / total_batch        if (time + 1) % train_step == 0:            accuracy = session.run(accuracy_op, feed_dict={input_x: mnist.validation.images})            summary_str = session.run(summary_op, feed_dict={input_x: mnist.validation.images})            summary_writer.add_summary(summary_str, session.run(index))            print("train times:", (time + 1),                        " avg_loss:", avg_loss,                        " accuracy:", accuracy)    # make to size 20 * 20 images form q(z) of distributed    imgs = np.empty((28 * 20, 28 * 20))    index_x = 0    for x in range(-10,10,1):        index_y = 0        for y in range(-10,10,1):            value = np.array([[float(x / 5.), float(y / 5.)]])            img = session.run(decoder_op, feed_dict={data_x:value})            imgs[index_x * 28:(index_x + 1) * 28, index_y * 28:(index_y + 1) * 28] = img.reshape((28, 28))            index_y += 1        index_x += 1    plt.imshow(imgs, cmap=plt.get_cmap('gray'))    plt.show()    session.close()

结果

结语

其实VAE整体上的公式理解并不算太複杂，而其实VAE最理想是求最大似然数，但因连续机率难以求出只好使用趋近方式，最后产生出的结果图片都是有些模糊的，可能是因为高斯分布的关係(在影像处理中的模糊也是其中用法)，接着介绍GAN，它可以使图片更加清晰的呈现，就让我们期待GAN会带来什么不同的惊喜吧。

补充：若对其中如何使用正态分布产生的图有兴趣可观看李弘毅老师的ML Lecture 18观看，里面也有另一个想法的公式推导。
另外这网站也有正态分布图片能观看，里面也是有不同的想法的公式推导。

参考文献

[1]cairohy(2018) 变分推断（一）：介绍变分推断 from: http://cairohy.github.io/2018/02/28/vi/VI-1/
[2]qixinbo(2018) 变分自编码器的原理和程序解析from: https://qixinbo.info/2018/07/24/vae/
[3]维基百科(2018) 相对熵 from: https://zh.wikipedia.org/wiki/%E7%9B%B8%E5%AF%B9%E7%86%B5
[4]冯超(2017) VAE(2)——基本思想 from: https://zhuanlan.zhihu.com/p/22464764
[5]极大似然估计详解(2017) 知行流浪 from: https://blog.csdn.net/zengxiantao1994/article/details/72787849
[6]维基百科(2018) 贝氏定理from: https://zh.wikipedia.org/wiki/%E8%B4%9D%E5%8F%B6%E6%96%AF%E5%AE%9A%E7%90%86