r-dl-cb-merge-1R 深度学习秘籍（二）第四章：使用自编码器的数据表示本章将介绍自编码器的无监督深度学习

R 深度学习秘籍（二）

原文：annas-archive.org/md5/fb6c5d96801d9512d19e799540d2d7ac

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：使用自编码器的数据表示

本章将介绍自编码器的无监督深度学习应用。我们将覆盖以下主题：

设置自编码器
数据归一化
设置正则化自编码器
微调自编码器的参数
设置堆叠自编码器
设置去噪自编码器
构建和比较随机编码器和解码器
从自编码器学习流形
评估稀疏分解

介绍

神经网络旨在找到输入X与输出y之间的非线性关系，即y=f(x)。自编码器是一种无监督神经网络，试图找到空间中特征之间的关系，使得h=f(x)，帮助我们学习输入空间之间的关系，可用于数据压缩、降维和特征学习。

自编码器由编码器和解码器组成。编码器帮助将输入x编码为潜在表示y，而解码器则将y转换回x。编码器和解码器具有相似的形式表示。

下面是单层自编码器的表示：

编码器将输入X编码为H，并通过隐藏层处理，而解码器则帮助从编码输出H恢复原始数据。矩阵W[e]和W[d]分别表示编码器和解码器层的权重。函数f是激活函数。

以下是自编码器的示意图：

以节点形式的约束使自编码器能够发现数据中的有趣结构。例如，在前面的图中，编码器的五个输入数据集必须经过三个节点的压缩才能得到编码值h。编码输出层的维度可以与输入/输出解码输出层相同、较低或较高。输入层节点数少于编码层的编码输出层称为欠完备表示，可以视为将数据压缩成低维表示。

编码输出层具有较多输入层的情况称为过完备表示，并在稀疏自编码器中作为正则化策略使用。自编码器的目标是找到y，捕捉数据变化的主要因素，这与**主成分分析（PCA）**类似，因此也可以用于压缩。

设置自编码器

存在多种不同架构的自编码器，这些架构通过使用不同的代价函数来捕捉数据表示。最基本的自编码器被称为香草自编码器（vanilla autoencoder）。它是一个包含两个层的神经网络，其中隐藏层的节点数与输入层和输出层相同，目标是最小化代价函数。常用的（但不限于）损失函数包括回归问题中的均方误差(MSE)和分类问题中的交叉熵。当前方法可以轻松扩展到多个层次，这种扩展被称为多层自编码器。

节点数量在自编码器中起着至关重要的作用。如果隐藏层中的节点数少于输入层的节点数，那么该自编码器被称为欠完备自编码器。隐藏层中节点数较多则表示过完备自编码器或稀疏自编码器。

稀疏自编码器旨在对隐藏层施加稀疏性。这种稀疏性可以通过在隐藏层引入比输入层更多的节点，或者通过在损失函数中引入惩罚来实现，从而使隐藏层的权重趋向于零。有些自编码器通过手动将节点的权重置为零来实现稀疏性，这些被称为K-稀疏自编码器。我们将在第一章《入门指南》中讨论的Occupancy数据集上设置自编码器。当前示例的隐藏层可以进行调整。

准备就绪

让我们使用Occupancy数据集来设置自编码器：

下载Occupancy数据集，如第一章《入门指南》中所述。
在 R 和 Python 中安装 TensorFlow

如何操作...

当前的Occupancy数据集，如第一章《入门指南》中所述，用于演示如何在 R 中使用 TensorFlow 设置自编码器：

设置 R TensorFlow 环境。
load_occupancy_data函数可以通过使用setwd设置正确的工作目录路径来加载数据：

# Function to load Occupancy data
load_occupancy_data<-function(train){
xFeatures = c("Temperature", "Humidity", "Light", "CO2",
  "HumidityRatio")
yFeatures = "Occupancy"
  if(train){
    occupancy_ds <-  as.matrix(read.csv("datatraining.txt",stringsAsFactors = T))
  } else
  {
    occupancy_ds <- as.matrix(read.csv("datatest.txt",stringsAsFactors = T))
  }
  occupancy_ds<-apply(occupancy_ds[, c(xFeatures, yFeatures)], 2, FUN=as.numeric) 
  return(occupancy_ds)
}

可以使用以下脚本将训练和测试的Occupancy数据集加载到 R 环境中：

occupancy_train <-load_occupancy_data(train=T)
occupancy_test <- load_occupancy_data(train = F)

数据归一化

数据归一化是机器学习中的一个关键步骤，用于将数据转换为相似的尺度。它也称为特征缩放，通常作为数据预处理的一部分进行。

正确的归一化对于神经网络至关重要，否则它会导致隐藏层的饱和，从而导致梯度为零，无法进行学习。

准备就绪

有多种方式可以执行归一化：

最小-最大标准化：最小-最大标准化保持原始分布，并将特征值缩放到*【0, 1】之间，其中0*为特征的最小值，1为最大值。标准化的过程如下：

这里，x' 是特征的归一化值。该方法对数据集中的异常值较为敏感。

十进制缩放：这种缩放形式用于存在不同十进制范围值的情况。例如，两个具有不同边界的特征可以通过如下方式使用十进制缩放将其带到相似的尺度：

x'=x/10^n

Z-得分：这种转换将值缩放到具有零均值和单位方差的正态分布。Z-得分计算公式为：

Z=(x-µ)/σ

这里，µ 是均值，σ 是特征的标准差。这些分布对具有高斯分布的数据集非常有效。

所有前述方法对异常值敏感；你可以探索其他更稳健的归一化方法，例如中位数绝对偏差 (MAD)、tanh 估计器和双重 sigmoid。

数据集分布的可视化

我们来看看职业数据的特征分布：

> ggpairs(occupancy_train$data[, occupancy_train$xFeatures])

图表显示特征具有线性相关性，并且分布是非正态的。可以通过使用 Shapiro-Wilk 测试进一步验证非正态性，使用 R 中的 shapiro.test 函数。我们将对职业数据使用最小-最大标准化。

如何操作...

执行以下操作以进行数据归一化：

minmax.normalize<-function(ds, scaler=NULL){
  if(is.null(scaler)){
    for(f in ds$xFeatures){
      scaler[[f]]$minval<-min(ds$data[,f])
      scaler[[f]]$maxval<-max(ds$data[,f])
      ds$data[,f]<-(ds$data[,f]-scaler[[f]]$minval)/(scaler[[f]]$maxval-scaler[[f]]$minval)
    }
    ds$scaler<-scaler
  } else
  {
    for(f in ds$xFeatures){
      ds$data[,f]<-(ds$data[,f]-scaler[[f]]$minval)/(scaler[[f]]$maxval-scaler[[f]]$minval)
    }
  }
  return(ds)
}

minmax.normalize 函数使用最小-最大标准化对数据进行归一化。当 scaler 变量为 NULL 时，它使用提供的数据集进行归一化，或者使用 scaler 的值进行归一化。归一化后的数据对图如下面的图所示：

此图显示了最小-最大归一化将值限制在 [0, 1] 范围内，并且没有改变特征之间的分布和相关性。

如何设置自编码器模型

下一步是设置自编码器模型。让我们使用 TensorFlow 设置一个基础的自编码器：

重置 graph 并启动 InteractiveSession：

# Reset the graph and set-up a interactive session
tf$reset_default_graph()
sess<-tf$InteractiveSession()

定义输入参数，其中 n 和 m 分别是样本数和特征数。为了构建网络，m 用来设置输入参数：

# Network Parameters
n_hidden_1 = 5 # 1st layer num features
n_input = length(xFeatures) # Number of input features
nRow<-nrow(occupancy_train)

当 n_hidden_1 很低时，自编码器会压缩数据，称为欠完备自编码器；而当 n_hidden_1 很大时，自编码器则是稀疏的，称为过完备自编码器。

定义包括输入张量和编码器、解码器层定义的图输入参数：

# Define input feature
x <- tf$constant(unlist(occupancy_train[, xFeatures]), shape=c(nRow, n_input), dtype=np$float32) 

# Define hidden and bias layer for encoder and decoders
hiddenLayerEncoder<-tf$Variable(tf$random_normal(shape(n_input, n_hidden_1)), dtype=np$float32)
biasEncoder <- tf$Variable(tf$zeros(shape(n_hidden_1)), dtype=np$float32)
hiddenLayerDecoder<-tf$Variable(tf$random_normal(shape(n_hidden_1, n_input)))
biasDecoder <- tf$Variable(tf$zeros(shape(n_input)))

前述脚本设计了一个单层编码器和解码器。

定义一个函数来评估响应：

auto_encoder<-function(x, hiddenLayerEncoder, biasEncoder){
  x_transform <- tf$nn$sigmoid(tf$add(tf$matmul(x, hiddenLayerEncoder), biasEncoder))
  x_transform
}

auto_encoder 函数接收节点偏置权重并计算输出。通过传递相应的权重，可以使用相同的函数来处理 encoder 和 decoder。

通过传递符号化的 TensorFlow 变量创建 encoder 和 decoder 对象：

encoder_obj = auto_encoder(x,hiddenLayerEncoder, biasEncoder)
y_pred = auto_encoder(encoder_obj, hiddenLayerDecoder, biasDecoder)

y_pred 是来自 decoder 的输出，它将 encoder 对象作为输入，包含节点和偏置权重：

Define loss function and optimizer module. 
learning_rate = 0.01
cost = tf$reduce_mean(tf$pow(x - y_pred, 2))
optimizer = tf$train$RMSPropOptimizer(learning_rate)$minimize(cost)

上述脚本将均方误差定义为成本函数，并使用 TensorFlow 中的RMSPropOptimizer，学习率为 0.1，来优化权重。上述模型的 TensorFlow 图如下所示：

运行优化

下一步是执行优化器优化。执行此过程的 TensorFlow 步骤包括两部分：

第一步是对图中定义的变量进行参数初始化。初始化通过调用 TensorFlow 中的global_variables_initializer函数来完成：

# Initializing the variables
init = tf$global_variables_initializer()
sess$run(init)

优化是基于优化和监控训练和测试性能进行的：

costconvergence<-NULL
for (step in 1:1000) {
  sess$run(optimizer)
  if (step %% 20==0){
    costconvergence<-rbind(costconvergence, c(step, sess$run(cost), sess$run(costt)))
    cat(step, "-", "Traing Cost ==>", sess$run(cost), "\n")
  }
}

可以观察训练和测试中的成本函数，以了解模型的收敛情况，如下图所示：

costconvergence<-data.frame(costconvergence)
colnames(costconvergence)<-c("iter", "train", "test")
plot(costconvergence[, "iter"], costconvergence[, "train"], type = "l", col="blue", xlab = "Iteration", ylab = "MSE")
lines(costconvergence[, "iter"], costconvergence[, "test"], col="red")
legend(500,0.25, c("Train","Test"), lty=c(1,1), lwd=c(2.5,2.5),col=c("blue","red"))

该图显示，模型的主要收敛发生在大约400次迭代时；然而，即使经过1,000次迭代后，收敛速度依然很慢。该模型在训练数据集和保留测试数据集中都很稳定。

设置正则化自编码器

正则化自编码器通过在cost函数中添加正则化参数来扩展标准自编码器。

准备工作

正则化自编码器是标准自编码器的扩展。设置将需要：

在 R 和 Python 中的 TensorFlow 安装。
标准自编码器的实现。

如何操作...

自编码器的代码设置可以通过将成本定义替换为以下几行，直接转换为正则化自编码器：

Lambda=0.01
cost = tf$reduce_mean(tf$pow(x - y_pred, 2))
Regularize_weights = tf$nn$l2_loss(weights)
cost = tf$reduce_mean(cost + lambda * Regularize_weights)

它是如何工作的...

如前所述，正则化自编码器通过在成本函数中添加正则化参数来扩展标准自编码器，如下所示：

这里，λ是正则化参数，i和j是节点索引，W表示自编码器的隐藏层权重。正则化自编码器的目的是确保更强大的编码，并偏好低权重h函数。这个概念进一步被用于开发收缩自编码器，它利用输入上雅可比矩阵的 Frobenius 范数，表示如下：

其中**J(x)**是雅可比矩阵，计算方法如下：

对于线性编码器，收缩编码器和正则化编码器收敛到 L2 权重衰减。正则化有助于使自编码器对输入的敏感度降低；然而，成本函数的最小化帮助模型捕捉变化，并保持对高密度流形的敏感性。这些自编码器也被称为收缩自编码器。

微调自编码器的参数

自编码器涉及一些需要调整的参数，这取决于我们所使用的自编码器类型。自编码器中的主要参数包括：

任何隐藏层中的节点数量
适用于深度自编码器的隐藏层数量
激活单元，例如 sigmoid、tanh、softmax 和 ReLU 激活函数
隐藏单元权重上的正则化参数或权重衰减项
在去噪自编码器中，信号损坏的比例
稀疏自编码器中的稀疏性参数，用于控制隐藏层中神经元的期望激活值
如果使用批量梯度下降学习，批次大小；如果使用随机梯度下降，学习率和动量参数
训练时使用的最大迭代次数
权重初始化
如果使用了 dropout，进行 dropout 正则化

这些超参数可以通过将问题设定为网格搜索问题来训练。然而，每个超参数组合都需要训练隐藏层的神经元权重，这导致随着层数和每层节点数的增加，计算复杂性也会增加。为了解决这些关键参数和训练问题，提出了堆叠自编码器概念，它逐层训练每一层，以获取预训练权重，然后使用获得的权重对模型进行微调。这种方法大大提高了训练性能，相比传统的训练方式。

设置堆叠自编码器

堆叠自编码器是一种训练深度网络的方法，包含多个层，使用贪婪算法逐层训练。下面的图示展示了堆叠自编码器的一个示例：

堆叠自编码器的示例

准备就绪

前面的图示展示了一个具有两层的堆叠自编码器。堆叠自编码器可以有n层，其中每一层是逐层训练的。例如，前一层的训练过程如下：

堆叠自编码器的训练

第 1 层的初步训练是通过在实际输入x[i]上进行训练得到的。第一步是优化编码器的We(1)层，以适应输出 X。上面的第二步是通过使用We(1)作为输入和输出，优化第二层的权重We(2)。一旦所有We(i)层（其中i=1, 2, ..., n为层数）都经过预训练，就可以通过将所有层连接在一起进行模型微调，正如前图中的第 3 步所示。该概念还可以应用于去噪训练多层网络，这被称为堆叠去噪自编码器。去噪自编码器中开发的代码可以轻松调整为开发堆叠去噪自编码器，它是堆叠自编码器的扩展。

本配方的要求如下：

需要安装 R。
SAENET包，可以通过命令install.packages("SAENET")从 Cran 下载安装该包。

如何操作...

在 R 中还有其他流行的库用于开发堆叠自编码器。让我们使用 R 中的SAENET包来搭建一个堆叠自编码器。SAENET是一个堆叠自编码器的实现，使用的是neuralnet包中的前馈神经网络（来自 CRAN）：

如果尚未安装，可以从 CRAN 仓库获取SAENET包：

install.packages("SAENET")

加载所有库依赖项：

require(SAENET)

使用load_occupancy_data加载训练和测试占用数据集：

occupancy_train <-load_occupancy_data(train=T)
occupancy_test <- load_occupancy_data(train = F)

使用minmax.normalize函数对数据集进行归一化：

# Normalize dataset
occupancy_train<-minmax.normalize(occupancy_train, scaler = NULL)
occupancy_test<-minmax.normalize(occupancy_test, scaler = occupancy_train$scaler)

堆叠自编码器模型可以使用SAENET.train训练函数从SAENET包中构建：

# Building Stacked Autoencoder
SAE_obj<-SAENET.train(X.train= subset(occupancy_train$data, select=-c(Occupancy)), n.nodes=c(4, 3, 2), unit.type ="tanh", lambda = 1e-5, beta = 1e-5, rho = 0.01, epsilon = 0.01, max.iterations=1000)

可以使用SAE_obj[[n]]$X.output命令提取最后一个节点的输出。

设置去噪自编码器

去噪自编码器是一种特殊类型的自编码器，专注于从输入数据集中提取稳健的特征。去噪自编码器与之前的模型类似，唯一的主要区别是数据在训练网络之前会被破坏。可以使用不同的破坏方法，例如遮罩，它会在数据中引入随机错误。

准备工作

让我们使用 CIFAR-10 图像数据来设置去噪数据集：

使用download_cifar_data函数下载 CIFAR-10 数据集（见第三章，卷积神经网络）
在 R 和 Python 中安装 TensorFlow

如何操作...

我们首先需要读取数据集。

读取数据集

使用第三章中解释的步骤加载CIFAR数据集，卷积神经网络。使用data_batch_1和data_batch_2文件进行训练，data_batch_5和test_batch文件分别用于验证和测试。数据可以通过flat_data函数进行展平：

train_data <- flat_data(x_listdata = images.rgb.train)
test_data <- flat_data(x_listdata = images.rgb.test)
valid_data <- flat_data(x_listdata = images.rgb.valid)

flat_data函数将数据集展平为NCOL = (高度 * 宽度 * 通道数)，因此数据集的维度是（图像数 × NCOL）。CIFAR中的图像大小为 32 x 32，包含三个 RGB 通道；因此，在数据展平后，我们得到 3,072 列：

> dim(train_data$images)
[1] 40000  3072

破坏数据以进行训练

设置去噪自编码器所需的下一个关键功能是数据破坏：

# Add noise using masking or salt & pepper noise method
add_noise<-function(data, frac=0.10, corr_type=c("masking", "saltPepper", "none")){
  if(length(corr_type)>1) corr_type<-corr_type[1] 

  # Assign a copy of data
  data_noise = data

  # Evaluate chaining parameters for autoencoder
  nROW<-nrow(data)
  nCOL<-ncol(data)
  nMask<-floor(frac*nCOL)

  if(corr_type=="masking"){
    for( i in 1:nROW){
      maskCol<-sample(nCOL, nMask)
      data_noise[i,maskCol,,]<-0
    }
  } else if(corr_type=="saltPepper"){
    minval<-min(data[,,1,])
    maxval<-max(data[,,1,])
    for( i in 1:nROW){
      maskCol<-sample(nCOL, nMask)
      randval<-runif(length(maskCol))
      ixmin<-randval<0.5
      ixmax<-randval>=0.5
      if(sum(ixmin)>0) data_noise[i,maskCol[ixmin],,]<-minval
      if(sum(ixmax)>0) data_noise[i,maskCol[ixmax],,]<-maxval
    }
  } else
  {
    data_noise<-data
  }
  return(data_noise)
}

可以使用以下脚本来破坏 CIFAR-10 数据：

# Corrupting input signal
xcorr<-add_noise(train_data$images, frac=0.10, corr_type="masking")

破坏后的示例图像如下：

上图使用了遮罩方法来添加噪声。此方法在随机图像位置添加零值，且定义了一个比例。另一种添加噪声的方法是使用椒盐噪声。此方法在图像中选择随机位置并进行替换，使用抛硬币原理为图像添加最小值或最大值。以下是使用椒盐方法进行数据破坏的示例：

数据破坏有助于自编码器学习更稳健的表示。

设置去噪自编码器

下一步是设置自编码器模型：

首先，重置图并开始一个交互式会话，如下所示：

# Reset the graph and set-up an interactive session
tf$reset_default_graph()
sess<-tf$InteractiveSession()

下一步是为输入信号和破坏信号定义两个占位符：

# Define Input as Placeholder variables
x = tf$placeholder(tf$float32, shape=shape(NULL, img_size_flat), name='x')
x_corrput<-tf$placeholder(tf$float32, shape=shape(NULL, img_size_flat), name='x_corrput')

x_corrupt 将作为自编码器的输入，而 x 是实际的图像，将用作输出。

设置去噪自编码器函数，如下代码所示：

# Setting-up denoising autoencoder
denoisingAutoencoder<-function(x, x_corrput, img_size_flat=3072, hidden_layer=c(1024, 512), out_img_size=256){

  # Building Encoder
  encoder = NULL
  n_input<-img_size_flat
  curentInput<-x_corrput
  layer<-c(hidden_layer, out_img_size)
  for(i in 1:length(layer)){
    n_output<-layer[i]
    W = tf$Variable(tf$random_uniform(shape(n_input, n_output), -1.0 / tf$sqrt(n_input), 1.0 / tf$sqrt(n_input)))
    b = tf$Variable(tf$zeros(shape(n_output)))
    encoder<-c(encoder, W)
    output = tf$nn$tanh(tf$matmul(curentInput, W) + b)
    curentInput = output
    n_input<-n_output
  }

  # latent representation
  z = curentInput
  encoder<-rev(encoder)
  layer_rev<-c(rev(hidden_layer), img_size_flat)

  # Build the decoder using the same weights
  decoder<-NULL
  for(i in 1:length(layer_rev)){
    n_output<-layer_rev[i]
    W = tf$transpose(encoder[[i]])
    b = tf$Variable(tf$zeros(shape(n_output)))
    output = tf$nn$tanh(tf$matmul(curentInput, W) + b)
    curentInput = output
  }

  # now have the reconstruction through the network
  y = curentInput

  # cost function measures pixel-wise difference
  cost = tf$sqrt(tf$reduce_mean(tf$square(y - x)))
  return(list("x"=x, "z"=z, "y"=y, "x_corrput"=x_corrput, "cost"=cost))
}

创建去噪对象：

# Create denoising AE object
dae_obj<-denoisingAutoencoder(x, x_corrput, img_size_flat=3072, hidden_layer=c(1024, 512), out_img_size=256)

设置代价函数：

# Learning set-up
learning_rate = 0.001
optimizer = tf$train$AdamOptimizer(learning_rate)$minimize(dae_obj$cost)

运行优化：

# We create a session to use the graph
sess$run(tf$global_variables_initializer())
for(i in 1:500){
  spls <- sample(1:dim(xcorr)[1],1000L)
  if (i %% 1 == 0) {
    x_corrput_ds<-add_noise(train_data$images[spls, ], frac = 0.3, corr_type = "masking")
    optimizer$run(feed_dict = dict(x=train_data$images[spls, ], x_corrput=x_corrput_ds))
    trainingCost<-dae_obj$cost$eval((feed_dict = dict(x=train_data$images[spls, ], x_corrput=x_corrput_ds)))
    cat("Training Cost - ", trainingCost, "\n")
  }
}

它是如何工作的...

自编码器继续学习关于特征的函数形式，以捕捉输入和输出之间的关系。计算机在 1,000 次迭代后如何可视化图像的示例如下所示：

在进行 1,000 次迭代后，计算机能够区分物体和环境的主要部分。随着我们继续运行算法来微调权重，计算机会继续学习更多关于物体本身的特征，如下图所示：

上图显示了模型仍在学习，但随着它开始学习关于物体的精细特征，学习率在迭代过程中逐渐变小，如下图所示。有时，模型开始上升而不是下降，这是由于批量梯度下降引起的：

使用去噪自编码器进行学习的示意图

构建和比较随机编码器和解码器

随机编码器属于生成建模领域，目标是学习给定数据 X 在转换到另一个高维空间后的联合概率 P(X)。例如，我们想学习图像并通过学习像素依赖关系和分布生成类似但不完全相同的图像。生成建模中一种流行的方法是变分自编码器（VAE），它通过对 h ~ P(h) 做出强假设，如高斯分布或伯努利分布，将深度学习与统计推理相结合。对于给定的权重 W，X 可以从分布中采样为 Pw(X|h)。下面的图示展示了 VAE 的架构：

VAE 的代价函数基于对数似然最大化。代价函数由重建误差和正则化误差项组成：

代价 = 重建误差 + 正则化误差

重建误差 是我们如何将结果与训练数据进行映射的准确度，而 正则化误差 对编码器和解码器形成的分布施加了惩罚。

准备工作

TensorFlow 需要在环境中安装并加载：

require(tensorflow)

需要加载依赖项：

require(imager)
require(caret)

MNIST 数据集需要被加载。数据集使用以下脚本进行归一化：

# Normalize Dataset
normalizeObj<-preProcess(trainData, method="range")
trainData<-predict(normalizeObj, trainData)
validData<-predict(normalizeObj, validData)

如何做...

MNIST 数据集被用来演示稀疏分解的概念。MNIST 数据集包含手写数字。它是从 tensorflow 数据集库下载的。数据集包含 28 x 28 像素的手写图像，包含 55,000 个训练样本，10,000 个测试样本和 5,000 个测试样本。可以通过以下脚本从 tensorflow 库下载数据集：

library(tensorflow)
datasets <- tf$contrib$learn$datasets
mnist <- datasets$mnist$read_data_sets("MNIST-data", one_hot = TRUE)

为了简化计算，MNIST 图像的大小从 28 x 28 像素减小为 16 x 16 像素，使用以下函数：

# Function to reduce image size
reduceImage<-function(actds, n.pixel.x=16, n.pixel.y=16){
  actImage<-matrix(actds, ncol=28, byrow=FALSE)
  img.col.mat <- imappend(list(as.cimg(actImage)),"c")
  thmb <- resize(img.col.mat, n.pixel.x, n.pixel.y)
  outputImage<-matrix(thmb[,,1,1], nrow = 1, byrow = F)
  return(outputImage)
}

以下脚本可以用来准备具有 16 x 16 像素图像的 MNIST 训练数据：

# Covert train data to 16 x 16  image
trainData<-t(apply(mnist$train$images, 1, FUN=reduceImage))
validData<-t(apply(mnist$test$images, 1, FUN=reduceImage))

plot_mnist 函数可以用来可视化选择的 MNIST 图像：

# Function to plot MNIST dataset
plot_mnist<-function(imageD, pixel.y=16){
  actImage<-matrix(imageD, ncol=pixel.y, byrow=FALSE)
  img.col.mat <- imappend(list(as.cimg(actImage)), "c")
  plot(img.col.mat)
}

设置 VAE 模型：

启动一个新的 TensorFlow 环境：

tf$reset_default_graph()
sess<-tf$InteractiveSession()

定义网络参数：

n_input=256
n.hidden.enc.1<-64

启动一个新的 TensorFlow 环境：

tf$reset_default_graph()
sess<-tf$InteractiveSession()

定义网络参数：

n_input=256
n.hidden.enc.1<-64

前述参数将形成如下的 VAE 网络：

定义模型初始化函数，定义每一层 encoder 和 decoder 的权重和偏差：

model_init<-function(n.hidden.enc.1, n.hidden.enc.2, 
                               n.hidden.dec.1,  n.hidden.dec.2, 
                               n_input, n_h)
{ weights<-NULL
 ############################
 # Set-up Encoder
 ############################
 # Initialize Layer 1 of encoder
 weights[["encoder_w"]][["h1"]]=tf$Variable(xavier_init(n_input,
 n.hidden.enc.1))
 weights[["encoder_w"]]
[["h2"]]=tf$Variable(xavier_init(n.hidden.enc.1, n.hidden.enc.2))
  weights[["encoder_w"]][["out_mean"]]=tf$Variable(xavier_init(n.hidden.enc.2, n_h))
  weights[["encoder_w"]][["out_log_sigma"]]=tf$Variable(xavier_init(n.hidden.enc.2, n_h))
  weights[["encoder_b"]][["b1"]]=tf$Variable(tf$zeros(shape(n.hidden.enc.1), dtype=tf$float32))
  weights[["encoder_b"]][["b2"]]=tf$Variable(tf$zeros(shape(n.hidden.enc.2), dtype=tf$float32))
  weights[["encoder_b"]][["out_mean"]]=tf$Variable(tf$zeros(shape(n_h), dtype=tf$float32))
  weights[["encoder_b"]][["out_log_sigma"]]=tf$Variable(tf$zeros(shape(n_h), dtype=tf$float32))

 ############################
 # Set-up Decoder
 ############################
 weights[['decoder_w']][["h1"]]=tf$Variable(xavier_init(n_h, n.hidden.dec.1))
 weights[['decoder_w']][["h2"]]=tf$Variable(xavier_init(n.hidden.dec.1, n.hidden.dec.2))
 weights[['decoder_w']][["out_mean"]]=tf$Variable(xavier_init(n.hidden.dec.2, n_input))
 weights[['decoder_w']][["out_log_sigma"]]=tf$Variable(xavier_init(n.hidden.dec.2, n_input))
 weights[['decoder_b']][["b1"]]=tf$Variable(tf$zeros(shape(n.hidden.dec.1), dtype=tf$float32))
 weights[['decoder_b']][["b2"]]=tf$Variable(tf$zeros(shape(n.hidden.dec.2), dtype=tf$float32))
 weights[['decoder_b']][["out_mean"]]=tf$Variable(tf$zeros(shape(n_input), dtype=tf$float32))
 weights[['decoder_b']][["out_log_sigma"]]=tf$Variable(tf$zeros(shape(n_input), dtype=tf$float32))
 return(weights)
}

model_init 函数返回 weights，它是一个二维列表。第一个维度表示权重的关联和类型。例如，它描述了 weights 变量是分配给编码器还是解码器，并且它是否存储节点的权重或偏差。model_init 中的 xavier_init 函数用于为模型训练分配初始权重：

# Xavier Initialization using Uniform distribution 
xavier_init<-function(n_inputs, n_outputs, constant=1){
  low = -constant*sqrt(6.0/(n_inputs + n_outputs)) 
  high = constant*sqrt(6.0/(n_inputs + n_outputs))
  return(tf$random_uniform(shape(n_inputs, n_outputs), minval=low, maxval=high, dtype=tf$float32))
}

设置编码器评估函数：

# Encoder update function
vae_encoder<-function(x, weights, biases){
  layer_1 = tf$nn$softplus(tf$add(tf$matmul(x, weights[['h1']]), biases[['b1']])) 
  layer_2 = tf$nn$softplus(tf$add(tf$matmul(layer_1, weights[['h2']]), biases[['b2']])) 
  z_mean = tf$add(tf$matmul(layer_2, weights[['out_mean']]), biases[['out_mean']])
  z_log_sigma_sq = tf$add(tf$matmul(layer_2, weights[['out_log_sigma']]), biases[['out_log_sigma']])
  return (list("z_mean"=z_mean, "z_log_sigma_sq"=z_log_sigma_sq))
}

vae_encoder 计算均值和方差，用于从隐藏层的权重和偏差中采样：

设置解码器评估函数：

# Decoder update function
vae_decoder<-function(z, weights, biases){
  layer1<-tf$nn$softplus(tf$add(tf$matmul(z, weights[["h1"]]), biases[["b1"]]))
  layer2<-tf$nn$softplus(tf$add(tf$matmul(layer1, weights[["h2"]]), biases[["b2"]]))
  x_reconstr_mean<-tf$nn$sigmoid(tf$add(tf$matmul(layer2, weights[['out_mean']]), biases[['out_mean']]))
  return(x_reconstr_mean)
}

vae_decoder 函数计算与采样层相关的均值和标准差，以及输出和平均输出：

设置重构估计的函数：

# Parameter evaluation
network_ParEval<-function(x, network_weights, n_h){

  distParameter<-vae_encoder(x, network_weights[["encoder_w"]], network_weights[["encoder_b"]])
  z_mean<-distParameter$z_mean
  z_log_sigma_sq <-distParameter$z_log_sigma_sq

  # Draw one sample z from Gaussian distribution
  eps = tf$random_normal(shape(BATCH, n_h), 0, 1, dtype=tf$float32)

  # z = mu + sigma*epsilon
  z = tf$add(z_mean, tf$multiply(tf$sqrt(tf$exp(z_log_sigma_sq)), eps))

  # Use generator to determine mean of
  # Bernoulli distribution of reconstructed input
  x_reconstr_mean <- vae_decoder(z, network_weights[["decoder_w"]], network_weights[["decoder_b"]])
  return(list("x_reconstr_mean"=x_reconstr_mean, "z_log_sigma_sq"=z_log_sigma_sq, "z_mean"=z_mean))
}

定义优化的成本函数：

# VAE cost function
vae_optimizer<-function(x, networkOutput){
  x_reconstr_mean<-networkOutput$x_reconstr_mean
  z_log_sigma_sq<-networkOutput$z_log_sigma_sq
  z_mean<-networkOutput$z_mean
  loss_reconstruction<--1*tf$reduce_sum(x*tf$log(1e-10 + x_reconstr_mean)+
                                       (1-x)*tf$log(1e-10 + 1 - x_reconstr_mean), reduction_indices=shape(1))
  loss_latent<--0.5*tf$reduce_sum(1+z_log_sigma_sq-tf$square(z_mean)-
                                    tf$exp(z_log_sigma_sq), reduction_indices=shape(1))
  cost = tf$reduce_mean(loss_reconstruction + loss_latent)
  return(cost)
}

设置模型进行训练：

# VAE Initialization
x = tf$placeholder(tf$float32, shape=shape(NULL, img_size_flat), name='x')
network_weights<-model_init(n.hidden.enc.1, n.hidden.enc.2, 
                            n.hidden.dec.1,  n.hidden.dec.2, 
                            n_input, n_h)
networkOutput<-network_ParEval(x, network_weights, n_h)
cost=vae_optimizer(x, networkOutput)
optimizer = tf$train$AdamOptimizer(lr)$minimize(cost)

运行优化：

sess$run(tf$global_variables_initializer())
for(i in 1:ITERATION){
  spls <- sample(1:dim(trainData)[1],BATCH)
  out<-optimizer$run(feed_dict = dict(x=trainData[spls,]))
  if (i %% 100 == 0){
  cat("Iteration - ", i, "Training Loss - ",  cost$eval(feed_dict = dict(x=trainData[spls,])), "\n")
  }
}

VAE 自编码器的输出：

可以使用以下脚本生成结果：

spls <- sample(1:dim(trainData)[1],BATCH)
networkOutput_run<-sess$run(networkOutput, feed_dict = dict(x=trainData[spls,]))

# Plot reconstructured Image
x_sample<-trainData[spls,]
NROW<-nrow(networkOutput_run$x_reconstr_mean)
n.plot<-5
par(mfrow = c(n.plot, 2), mar = c(0.2, 0.2, 0.2, 0.2), oma = c(3, 3, 3, 3))
pltImages<-sample(1:NROW,n.plot)
for(i in pltImages){
  plot_mnist(x_sample[i,])
  plot_mnist(networkOutput_run$x_reconstr_mean[i,])
}

以下图所示是从前述 VAE 自编码器经过 20,000 次迭代后的结果：

此外，由于 VAE 是生成模型，输出结果不是输入的精确复制，而是会随着运行的不同而有所变化，因为从估计的分布中提取了一个代表性样本。

从自编码器学习流形：

流形学习是机器学习中的一种方法，它假设数据位于一个维度远低于原始数据的流形上。这些流形可以是线性或非线性的。因此，该方法尝试将数据从高维空间投影到低维空间。例如，主成分分析（PCA）是线性流形学习的一个例子，而自编码器则是一个非线性降维（NDR）方法，具有学习低维空间中非线性流形的能力。线性与非线性流形学习的比较见下图：

如图**a)所示，数据位于一个线性流形上，而在图b)**中，数据则位于二阶非线性流形上。

如何操作...

让我们取堆叠自编码器部分的输出，分析数据在转移到不同维度时，流形的表现。

设置主成分分析

在进入非线性流形之前，让我们分析占用数据上的主成分分析：

# Setting-up principal component analysis 
pca_obj <- prcomp(occupancy_train$data,
                 center = TRUE,
                 scale. = TRUE)
                 scale. = TRUE)

上述函数将数据转化为六个正交方向，这些方向是特征的线性组合。每个维度所解释的方差可以通过以下脚本查看：

plot(pca_obj, type = "l")

上述命令将绘制主成分的方差，如下图所示：

对于占用数据集，前两个主成分捕捉了大部分变异，当绘制主成分时，显示出占用的正负类之间的分离，如下图所示：

前两个主成分的输出

让我们可视化自编码器学习到的低维流形。我们只使用一个维度来可视化结果，如下所示：

SAE_obj<-SAENET.train(X.train= subset(occupancy_train$data, select=-c(Occupancy)), n.nodes=c(4, 3, 1), unit.type ="tanh", lambda = 1e-5, beta = 1e-5, rho = 0.01, epsilon = 0.01, max.iterations=1000)

上述脚本的编码器架构如下所示：

堆叠自编码器中一个潜在节点的隐藏层输出如下所示：

上述图表显示，占用在潜在变量的峰值处为真。然而，峰值出现在不同的值上。让我们增加潜在变量 2，这是通过 PCA 捕获的。可以使用以下脚本开发模型并绘制数据：

SAE_obj<-SAENET.train(X.train= subset(occupancy_train$data, select=-c(Occupancy)), n.nodes=c(4, 3, 2), unit.type ="tanh", lambda = 1e-5, beta = 1e-5, rho = 0.01, epsilon = 0.01, max.iterations=1000)  

# plotting encoder values
plot(SAE_obj[[3]]$X.output[,1], SAE_obj[[3]]$X.output[,2], col="blue", xlab = "Node 1 of layer 3", ylab = "Node 2 of layer 3")
ix<-occupancy_train$data[,6]==1  
points(SAE_obj[[3]]$X.output[ix,1], SAE_obj[[3]]$X.output[ix,2], col="red")

带有两层编码的值在下图中显示：

评估稀疏分解

稀疏自编码器也被称为过完备表示，并且在隐藏层中具有更多的节点。稀疏自编码器通常使用稀疏性参数（正则化）执行，该参数充当约束，限制节点的激活。稀疏性也可以被看作是由于稀疏性约束而导致的节点丢弃。稀疏自编码器的损失函数由重建误差、用于控制权重衰减的正则化项以及用于稀疏性约束的 KL 散度组成。以下表示很好地说明了我们所讲的内容：

准备工作

数据集已加载并设置好。
使用以下脚本安装并加载autoencoder包：

install.packages("autoencoder")
require(autoencoder)

如何操作...

TensorFlow 的标准自编码器代码可以通过更新成本函数轻松扩展到稀疏自编码器模块。本节将介绍 R 的自编码器包，该包内置了运行稀疏自编码器的功能：

### Setting-up parameter
nl<-3 
N.hidden<-100 
unit.type<-"logistic" 
lambda<-0.001 
rho<-0.01 
beta<-6 
max.iterations<-2000 
epsilon<-0.001 

### Running sparse autoencoder
spe_ae_obj <- autoencode(X.train=trainData,  X.test = validData, nl=nl, N.hidden=N.hidden, unit.type=unit.type,lambda=lambda,beta=beta,              epsilon=epsilon,rho=rho,max.iterations=max.iterations, rescale.flag = T)

autoencode函数中的主要参数如下：

nl：这是包括输入层和输出层在内的层数（默认值为三层）。
N.hidden：这是每个隐藏层中神经元数量的向量。
unit.type：这是要使用的激活函数类型。
lambda：这是正则化参数。
rho：这是稀疏性参数。
beta：这是稀疏项的惩罚。
max.iterations：这是最大迭代次数。
epsilon：这是权重初始化的参数。权重是使用高斯分布 ~N(0, epsilon2) 初始化的。

它是如何工作的...

下图展示了稀疏自编码器捕获的来自MNIST的数字形状和方向：

通过稀疏自编码器生成的滤波器得到数字结果

稀疏自编码器学习到的滤波器可以使用来自自编码器包的visualize.hidden.units函数进行可视化。该包绘制了最终层的权重与输出之间的关系。在当前场景中，100 是隐藏层中的神经元数量，256 是输出层中的节点数量。

第五章：深度学习中的生成模型

在本章中，我们将讨论以下主题：

比较主成分分析与限制玻尔兹曼机
为伯努利分布输入设置限制玻尔兹曼机
训练限制玻尔兹曼机
RBM 的反向或重建阶段
理解重建的对比散度
初始化并启动一个新的 TensorFlow 会话
评估 RBM 的输出
为协同过滤设置限制玻尔兹曼机
执行 RBM 训练的完整运行
设置深度信念网络
实现前馈反向传播神经网络
设置深度限制玻尔兹曼机

比较主成分分析与限制玻尔兹曼机

在本节中，你将学习两种广泛推荐的降维技术——主成分分析（PCA）和限制玻尔兹曼机（RBM）。考虑在n维空间中的一个向量v。降维技术本质上将向量v转换为一个相对较小（有时是相等）的m维向量v'（m<n）。这种转换可以是线性或非线性的。

PCA 对特征进行线性变换，从而生成正交调整的组件，这些组件之后会根据它们在方差捕捉中的相对重要性进行排序。这些m个组件可以视为新的输入特征，并可以如下定义：

向量v' =

这里，w 和 c 分别对应于权重（加载）和转换后的组件。

与 PCA 不同，RBM（或 DBN/自编码器）通过可见单元和隐藏单元之间的连接执行非线性变换，正如在第四章 使用自编码器的数据表示中所描述的那样。非线性有助于更好地理解与潜在变量之间的关系。除了信息捕获外，它们还倾向于去除噪声。RBM 通常基于随机分布（无论是伯努利分布还是高斯分布）。

执行大量的吉布斯采样以学习和优化可见层与隐藏层之间的连接权重。优化过程分为两个阶段：前向阶段，其中使用给定的可见层对隐藏层进行采样；反向阶段，其中使用给定的隐藏层对可见层进行重新采样。该优化旨在最小化重建误差。

以下图像表示一个限制玻尔兹曼机：

准备工作

对于这个食谱，你将需要 R（rbm和ggplot2包）和 MNIST 数据集。MNIST 数据集可以从 TensorFlow 数据集库中下载。该数据集包含 28 x 28 像素的手写图像。它有 55,000 个训练样本和 10,000 个测试样本。可以通过以下脚本从tensorflow库下载：

library(tensorflow) 
datasets <- tf$contrib$learn$datasets 
mnist <- datasets$mnist$read_data_sets("MNIST-data", one_hot = TRUE)

如何操作...

提取训练数据集（trainX包含所有 784 个独立变量，trainY包含相应的 10 个二元输出）：

trainX <- mnist$train$images 
trainY <- mnist$train$labels

对trainX数据执行 PCA：

PCA_model <- prcomp(trainX, retx=TRUE)

对trainX数据运行 RBM：

RBM_model <- rbm(trainX, retx=TRUE, max_epoch=500,num_hidden =900)

使用生成的模型对训练数据进行预测。对于 RBM 模型，生成概率：

PCA_pred_train <- predict(PCA_model) 
RBM_pred_train <- predict(RBM_model,type='probs')

将结果转换为数据框：

PCA_pred_train <- as.data.frame(PCA_pred_train)
 class="MsoSubtleEmphasis">RBM_pred_train <- as.data.frame(as.matrix(RBM_pred_train))

将 10 类二元trainY数据框转换为数值向量：

    trainY_num<- as.numeric(stringi::stri_sub(colnames(as.data.frame(trainY))[max.col(as.data.frame(trainY),ties.method="first")],2))

绘制使用 PCA 生成的组件图。这里，x轴表示组件 1，y轴表示组件 2。以下图片展示了 PCA 模型的结果：

ggplot(PCA_pred_train, aes(PC1, PC2))+
  geom_point(aes(colour = trainY))+ 
  theme_bw()+labs()+ 
  theme(plot.title = element_text(hjust = 0.5))

绘制使用 PCA 生成的隐藏层。这里，x轴表示隐藏 1，y轴表示隐藏 2。以下图片展示了 RBM 模型的结果：

ggplot(RBM_pred_train, aes(Hidden_2, Hidden_3))+
  geom_point(aes(colour = trainY))+ 
  theme_bw()+labs()+ 
  theme(plot.title = element_text(hjust = 0.5))

以下代码和图片展示了主成分所解释的累积方差：

var_explain <- as.data.frame(PCA_model$sdev²/sum(PCA_model$sdev²)) 
var_explain <- cbind(c(1:784),var_explain,cumsum(var_explain[,1])) 
colnames(var_explain) <- c("PcompNo.","Ind_Variance","Cum_Variance") 
plot(var_explain$PcompNo.,var_explain$Cum_Variance, xlim = c(0,100),type='b',pch=16,xlab = "# of Principal Components",ylab = "Cumulative Variance",main = 'PCA - Explained variance')

以下代码和图片展示了在使用多个训练周期生成 RBM 时，重建训练误差的下降：

plot(RBM_model,xlab = "# of epoch iterations",ylab = "Reconstruction error",main = 'RBM - Reconstruction Error')

为伯努利分布输入设置限制玻尔兹曼机

在本节中，我们将为伯努利分布的输入数据设置限制玻尔兹曼机，其中每个属性的值范围从 0 到 1（相当于一个概率分布）。本配方中使用的数据集（MNIST）具有满足伯努利分布的输入数据。

限制玻尔兹曼机由两层组成：一个可见层和一个隐藏层。可见层是输入层，节点数量等于输入属性的数量。在我们的案例中，MNIST 数据集中的每个图像由 784 个像素（28 x 28 大小）定义。因此，我们的可见层将包含 784 个节点。

另一方面，隐藏层通常是由用户定义的。隐藏层具有一组二值激活的节点，每个节点与所有其他可见节点有一定的连接概率。在我们的案例中，隐藏层将包含 900 个节点。作为初步步骤，所有可见层的节点与所有隐藏层的节点是双向连接的。

每个连接都由一个权重定义，因此定义了一个权重矩阵，其中行代表输入节点的数量，列代表隐藏节点的数量。在我们的案例中，权重矩阵（w）将是一个 784 x 900 的张量。

除了权重外，每个层中的所有节点还由偏置节点辅助。可见层的偏置节点将与所有可见节点（即 784 个节点）连接，表示为vb，而隐藏层的偏置节点将与所有隐藏节点（即 900 个节点）连接，表示为vh。

记住，RBM 的一个要点是每层内部的节点之间没有连接。换句话说，连接是跨层的，而不是层内的。

以下图像表示了包含可见层、隐藏层和连接的 RBM：

准备工作

本节提供了设置 RBM 的要求。

在 R 中安装并设置 TensorFlow
mnist数据被下载并加载以设置 RBM

如何操作...

本节提供了使用 TensorFlow 设置 RBM 的可见层和隐藏层的步骤：

启动一个新的交互式 TensorFlow 会话：

# Reset the graph 
tf$reset_default_graph() 
# Starting session as interactive session 
sess <- tf$InteractiveSession()

定义模型参数。num_input参数定义可见层节点的数量，num_hidden定义隐藏层节点的数量：

num_input<-784L 
num_hidden<-900L

为权重矩阵创建占位符变量：

W <- tf$placeholder(tf$float32, shape = shape(num_input, num_hidden))

为可见和隐藏偏置创建占位符变量：

vb <- tf$placeholder(tf$float32, shape = shape(num_input)) 
hb <- tf$placeholder(tf$float32, shape = shape(num_hidden))

训练限制玻尔兹曼机

每次训练 RBM 都会经历两个阶段：正向阶段和反向阶段（或重建阶段）。通过进行正向和反向阶段的多次迭代，来精细调节可见单元的重建。

正向阶段训练：在正向阶段，输入数据从可见层传递到隐藏层，所有计算发生在隐藏层的节点中。计算本质上是对每个从可见层到隐藏层连接的随机决策。在隐藏层中，输入数据（X）与权重矩阵（W）相乘，并加上一个隐藏偏置向量（hb）。

得到的向量（大小等于隐藏层节点的数量）会通过 sigmoid 函数，确定每个隐藏节点的输出（或激活状态）。在我们的案例中，每个输入数字会生成一个包含 900 个概率的张量向量，由于我们有 55,000 个输入数字，所以我们会得到一个大小为 55,000 x 900 的激活矩阵。利用隐藏层的概率分布矩阵，我们可以生成激活向量的样本，之后可用于估计负向阶段的梯度。

准备工作

本节提供了设置 RBM 的要求。

在 R 中安装并设置 TensorFlow
mnist数据被下载并加载以设置 RBM
RBM 模型的设置按照为伯努利分布输入设置限制玻尔兹曼机的步骤进行。

一个采样的示例

假设一个常量向量s1，它等于一个概率的张量向量。然后，使用常量向量s1的分布，创建一个新的随机均匀分布样本s2。接着计算它们的差异，并应用一个修正的线性激活函数。

如何操作...

本节提供了使用 TensorFlow 运行 RBM 模型的脚本设置步骤：

X = tf$placeholder(tf$float32, shape=shape(NULL, num_input)) 
prob_h0= tf$nn$sigmoid(tf$matmul(X, W) + hb) 
h0 = tf$nn$relu(tf$sign(prob_h0 - tf$random_uniform(tf$shape(prob_h0))))

使用以下代码执行在 TensorFlow 中创建的图：

sess$run(tf$global_variables_initializer()) 
s1 <- tf$constant(value = c(0.1,0.4,0.7,0.9)) 
cat(sess$run(s1)) 
s2=sess$run(tf$random_uniform(tf$shape(s1))) 
cat(s2) 
cat(sess$run(s1-s2)) 
cat(sess$run(tf$sign(s1 - s2))) 
cat(sess$run(tf$nn$relu(tf$sign(s1 - s2))))

RBM 的反向或重建阶段

在重建阶段，来自隐藏层的数据被传递回可见层。隐藏层的概率向量h0与权重矩阵W的转置相乘，并加上一个可见层偏置vb，然后通过 Sigmoid 函数生成重建的输入向量prob_v1。

通过使用重建的输入向量创建一个样本输入向量，该向量随后与权重矩阵W相乘，并加上隐藏偏置向量hb，生成更新后的隐藏概率向量h1。

这也叫做吉布斯采样。在某些情况下，样本输入向量不会生成，而是直接使用重建的输入向量prob_v1来更新隐藏层

准备开始

本节提供了使用输入概率向量进行图像重建的要求。

mnist数据已加载到环境中
RBM 模型是通过训练限制玻尔兹曼机的配方进行训练的

如何做...

本节介绍了执行反向重建和评估的步骤：

反向图像重建可以使用输入概率向量通过以下脚本进行：

prob_v1 = tf$nn$sigmoid(tf$matmul(h0, tf$transpose(W)) + vb)
v1 = tf$nn$relu(tf$sign(prob_v1 - tf$random_uniform(tf$shape(prob_v1))))
h1 = tf$nn$sigmoid(tf$matmul(v1, W) + hb)

评估可以使用一个定义的度量标准进行，如均方误差（MSE），它是在实际输入数据（X）和重建的输入数据（v1）之间计算的。MSE 在每个周期后计算，关键目标是最小化 MSE：

err = tf$reduce_mean(tf$square(X - v1))

理解重建的对比散度

作为初步设置，目标函数可以定义为最小化重建可见向量v的平均负对数似然，其中*P(v)*表示生成概率的向量：

准备开始

本节提供了使用输入概率向量进行图像重建的要求。

mnist数据已加载到环境中
图像是使用反向或重建阶段的配方重建的

如何做...

当前配方介绍了使用对比散度（CD）技术加速采样过程的步骤：

通过将输入向量X与来自给定概率分布prob_h0的隐藏向量h0样本相乘（外积），计算正权重梯度：

w_pos_grad = tf$matmul(tf$transpose(X), h0)

通过将重建的输入数据样本v1与更新的隐藏激活向量h1进行外积，计算负权重梯度：

w_neg_grad = tf$matmul(tf$transpose(v1), h1)

然后，通过从正梯度中减去负梯度并除以输入数据的大小来计算CD矩阵：

CD = (w_pos_grad - w_neg_grad) / tf$to_float(tf$shape(X)[0])

然后，使用学习率（alpha）和 CD 矩阵，将权重矩阵W更新为update_W：

update_w = W + alpha * CD

此外，更新可见和隐藏偏置向量：

update_vb = vb + alpha * tf$reduce_mean(X - v1) 
update_hb = hb + alpha * tf$reduce_mean(h0 - h1)

它是如何工作的...

目标函数可以通过随机梯度下降法最小化，间接修改（并优化）权重矩阵。整个梯度可以基于概率密度进一步分为两种形式：正梯度和负梯度。正梯度主要依赖于输入数据，而负梯度仅依赖于生成的模型。

在正梯度中，重构训练数据的概率增加，而在负梯度中，由模型随机生成的均匀样本的概率减小。

CD 技术用于优化负相。使用 CD 技术时，在每次重构迭代中都会调整权重矩阵。新的权重矩阵通过以下公式生成。学习率定义为alpha，在我们的情况下：

初始化并启动新的 TensorFlow 会话

计算误差度量（如均方误差 MSE）的一个重要部分是初始化和启动新的 TensorFlow 会话。以下是我们进行操作的步骤。

准备就绪

本节提供了启动新的 TensorFlow 会话所需的要求，用于计算误差度量。

mnist 数据已加载到环境中
RBM 的 TensorFlow 图已加载

如何操作...

本节提供了使用 RBM 重构优化误差的步骤：

初始化当前和前一个偏置向量及权重矩阵：

cur_w = tf$Variable(tf$zeros(shape = shape(num_input, num_hidden), dtype=tf$float32)) 
cur_vb = tf$Variable(tf$zeros(shape = shape(num_input), dtype=tf$float32)) 
cur_hb = tf$Variable(tf$zeros(shape = shape(num_hidden), dtype=tf$float32)) 
prv_w = tf$Variable(tf$random_normal(shape=shape(num_input, num_hidden), stddev=0.01, dtype=tf$float32)) 
prv_vb = tf$Variable(tf$zeros(shape = shape(num_input), dtype=tf$float32)) 
prv_hb = tf$Variable(tf$zeros(shape = shape(num_hidden), dtype=tf$float32))

启动一个新的 TensorFlow 会话：

sess$run(tf$global_variables_initializer())

使用完整的输入数据（trainX）执行第一次运行，并获得第一组权重矩阵和偏置向量：

output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(X=trainX, 
                                                                          W = prv_w$eval(), 
                                                                          vb = prv_vb$eval(), 
                                                                          hb = prv_hb$eval())) 
prv_w <- output[[1]]
prv_vb <-output[[2]] 
prv_hb <-output[[3]]

我们来看一下第一次运行的误差：

sess$run(err, feed_dict=dict(X= trainX, W= prv_w, vb= prv_vb, hb= prv_hb))

使用以下脚本可以训练 RBM 的完整模型：

epochs=15 
errors <- list() 
weights <- list() 
u=1 
for(ep in 1:epochs){ 
  for(i in seq(0,(dim(trainX)[1]-100),100)){ 
    batchX <- trainX[(i+1):(i+100),] 
    output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(X=batchX, 
                                                                              W = prv_w, 
                                                                              vb = prv_vb, 
                                                                              hb = prv_hb)) 
    prv_w <- output[[1]]  
    prv_vb <- output[[2]] 
    prv_hb <-  output[[3]] 
    if(i%%10000 == 0){ 
      errors[[u]] <- sess$run(err, feed_dict=dict(X= trainX, W= prv_w, vb= prv_vb, hb= prv_hb)) 
      weights[[u]] <- output[[1]] 
      u <- u+1 
     cat(i , " : ") 
    } 
  } 
  cat("epoch :", ep, " : reconstruction error : ", errors[length(errors)][[1]],"\n") 
}

使用均方误差绘制重构：

error_vec <- unlist(errors)
plot(error_vec,xlab="# of batches",ylab="mean squared reconstruction error",main="RBM-Reconstruction MSE plot")

它是如何工作的...

在这里，我们将运行 15 个周期（或迭代），每个周期中会执行批量（大小=100）的优化。在每个批次中，计算 CD 并相应更新权重和偏置。为了跟踪优化过程，每处理完 10,000 行数据后都会计算一次 MSE。

下图展示了计算 90 个批次的均方重构误差的下降趋势：

评估来自 RBM 的输出

在这里，我们将绘制最终层的权重与输出（重构输入数据）之间的关系。在当前的场景下，900 是隐藏层的节点数，784 是输出（重构）层的节点数。

在下图中，可以看到隐藏层的前 400 个节点：

在这里，每个图块表示一个隐藏节点与所有可见层节点之间的连接向量。在每个图块中，黑色区域表示负权重（权重 < 0），白色区域表示正权重（权重 > 1），灰色区域表示没有连接（权重 = 0）。正值越高，隐藏节点的激活概率越大，反之亦然。这些激活有助于确定给定隐藏节点正在确定输入图像的哪一部分。

准备开始

本节提供了运行评估教程所需的要求：

mnist数据已加载到环境中
使用 TensorFlow 执行 RBM 模型，并获得最佳权重

如何操作...

本教程涵盖了从 RBM 中获得的权重评估步骤：

运行以下代码生成 400 个隐藏节点的图像：

uw = t(weights[[length(weights)]])   # Extract the most recent weight matrix 
numXpatches = 20    # Number of images in X-axis (user input) 
numYpatches=20      # Number of images in Y-axis (user input) 
pixels <- list() 
op <- par(no.readonly = TRUE) 
par(mfrow = c(numXpatches,numYpatches), mar = c(0.2, 0.2, 0.2, 0.2), oma = c(3, 3, 3, 3)) 
for (i in 1:(numXpatches*numYpatches)) { 
  denom <- sqrt(sum(uw[i, ]²)) 
  pixels[[i]] <- matrix(uw[i, ]/denom, nrow = numYpatches, ncol = numXpatches) 
  image(pixels[[i]], axes = F, col = gray((0:32)/32)) 
} 
par(op)

从训练数据中选择四个实际输入的数字样本：

sample_image <- trainX[1:4,]

然后，使用以下代码可视化这些样本数字：

mw=melt(sample_image) 
mw$X3=floor((mw$X2-1)/28)+1 
mw$X2=(mw$X2-1)%%28 + 1; 
mw$X3=29-mw$X3 
ggplot(data=mw)+geom_tile(aes(X2,X3,fill=value))+facet_wrap(~X1,nrow=2)+ 
  scale_fill_continuous(low='black',high='white')+coord_fixed(ratio=1)+ 
  labs(x=NULL,y=NULL,)+ 
  theme(legend.position="none")+ 
  theme(plot.title = element_text(hjust = 0.5))

现在，使用最终获得的权重和偏置重构这四个样本图像：

hh0 = tf$nn$sigmoid(tf$matmul(X, W) + hb) 
vv1 = tf$nn$sigmoid(tf$matmul(hh0, tf$transpose(W)) + vb) 
feed = sess$run(hh0, feed_dict=dict( X= sample_image, W= prv_w, hb= prv_hb)) 
rec = sess$run(vv1, feed_dict=dict( hh0= feed, W= prv_w, vb= prv_vb))

然后，使用以下代码可视化重构后的样本数字：

mw=melt(rec) 
mw$X3=floor((mw$X2-1)/28)+1 
mw$X2=(mw$X2-1)%%28 + 1 
mw$X3=29-mw$X3 
ggplot(data=mw)+geom_tile(aes(X2,X3,fill=value))+facet_wrap(~X1,nrow=2)+ 
  scale_fill_continuous(low='black',high='white')+coord_fixed(ratio=1)+ 
  labs(x=NULL,y=NULL,)+ 
  theme(legend.position="none")+ 
  theme(plot.title = element_text(hjust = 0.5))

它是如何工作的...

下图展示了四个样本数字的原始图像：

重构的图像似乎去除了噪声，尤其是在数字3和6的情况下。

下图展示了同四个数字的重构图像：

为协同过滤设置限制玻尔兹曼机

在本教程中，你将学习如何使用 RBM 构建一个基于协同过滤的推荐系统。这里，对于每个用户，RBM 会根据他们过去对各种项目的评分行为，识别出相似的用户，然后尝试推荐下一个最佳项目。

准备开始

在本教程中，我们将使用 Grouplens 研究机构提供的 movielens 数据集。数据集（movies.dat 和 ratings.dat）可以通过以下链接下载。Movies.dat包含 3,883 部电影的信息，Ratings.dat包含 1,000,209 个用户对这些电影的评分。评分范围从 1 到 5，5 为最高分。

files.grouplens.org/datasets/movielens/ml-1m.zip

如何操作...

本教程涵盖了使用 RBM 设置协同过滤的步骤。

在 R 中读取movies.dat数据集：

txt <- readLines("movies.dat", encoding = "latin1") 
txt_split <- lapply(strsplit(txt, "::"), function(x) as.data.frame(t(x), stringsAsFactors=FALSE)) 
movies_df <- do.call(rbind, txt_split) 
names(movies_df) <- c("MovieID", "Title", "Genres") 
movies_df$MovieID <- as.numeric(movies_df$MovieID)

向电影数据集添加一个新列（id_order），因为当前的 ID 列（UserID）不能用来索引电影，因为它的值范围是从 1 到 3,952：

movies_df$id_order <- 1:nrow(movies_df)

在 R 中读取ratings.dat数据集：

ratings_df <- read.table("ratings.dat", sep=":",header=FALSE,stringsAsFactors = F) 
ratings_df <- ratings_df[,c(1,3,5,7)] 
colnames(ratings_df) <- c("UserID","MovieID","Rating","Timestamp")

使用all=FALSE合并电影和评分数据集：

merged_df <- merge(movies_df, ratings_df, by="MovieID",all=FALSE)

删除不必要的列：

merged_df[,c("Timestamp","Title","Genres")] <- NULL

将评分转换为百分比：

merged_df$rating_per <- merged_df$Rating/5

生成一个包含 1,000 个用户对所有电影评分的矩阵：

num_of_users <- 1000 
num_of_movies <- length(unique(movies_df$MovieID)) 
trX <- matrix(0,nrow=num_of_users,ncol=num_of_movies) 
for(i in 1:num_of_users){ 
  merged_df_user <- merged_df[merged_df$UserID %in% i,] 
  trX[i,merged_df_user$id_order] <- merged_df_user$rating_per 
}

查看trX训练数据集的分布。它似乎遵循伯努利分布（值的范围在 0 到 1 之间）：

summary(trX[1,]); summary(trX[2,]); summary(trX[3,])

定义输入模型参数：

num_hidden = 20 
num_input = nrow(movies_df)

启动一个新的 TensorFlow 会话：

sess$run(tf$global_variables_initializer()) 
output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(v0=trX, 
                                                                          W = prv_w$eval(), 
                                                                          vb = prv_vb$eval(), 
                                                                          hb = prv_hb$eval())) 
prv_w <- output[[1]]  
prv_vb <- output[[2]] 
prv_hb <-  output[[3]] 
sess$run(err_sum, feed_dict=dict(v0=trX, W= prv_w, vb= prv_vb, hb= prv_hb))

使用 500 个 epoch 迭代和批次大小 100 训练 RBM：

epochs= 500 
errors <- list() 
weights <- list() 

for(ep in 1:epochs){ 
  for(i in seq(0,(dim(trX)[1]-100),100)){ 
    batchX <- trX[(i+1):(i+100),] 
    output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(v0=batchX, 
                                                                              W = prv_w, 
                                                                              vb = prv_vb, 
                                                                              hb = prv_hb)) 
    prv_w <- output[[1]]  
    prv_vb <- output[[2]] 
    prv_hb <-  output[[3]] 
    if(i%%1000 == 0){ 
      errors <- c(errors,sess$run(err_sum, feed_dict=dict(v0=batchX, W= prv_w, vb= prv_vb, hb= prv_hb))) 
      weights <- c(weights,output[[1]]) 
      cat(i , " : ") 
    } 
  } 
  cat("epoch :", ep, " : reconstruction error : ", errors[length(errors)][[1]],"\n") 
}

绘制重建均方误差：

error_vec <- unlist(errors) 
plot(error_vec,xlab="# of batches",ylab="mean squared reconstruction error",main="RBM-Reconstruction MSE plot")

执行完整的 RBM 训练过程

使用前面配方中提到的相同 RBM 设置，使用 20 个隐藏节点对用户评分数据集（trX）进行训练。为了跟踪优化过程，每训练完 1,000 行就会计算一次 MSE。以下图片展示了计算得出的 500 批次（即 epochs）均方重建误差的下降趋势：

查看 RBM 推荐：现在让我们来看一下基于 RBM 的协同过滤为特定用户 ID 生成的推荐。在这里，我们将查看该用户 ID 的顶级评分流派和顶级推荐流派，并列出前 10 名电影推荐。

以下图片展示了顶级评分的流派列表：

以下图片展示了顶级推荐流派列表：

准备就绪

本节提供了协同过滤输出评估的要求：

已安装并设置 TensorFlow for R
movies.dat和ratings.dat数据集已加载到环境中
配方为协同过滤设置限制玻尔兹曼机已执行

如何执行...

本配方包含评估 RBM 协同过滤输出的步骤：

选择一个用户的评分：

inputUser = as.matrix(t(trX[75,]))
names(inputUser) <- movies_df$id_order

移除用户未评分的电影（假设这些电影还未观看）：

inputUser <- inputUser[inputUser>0]

绘制用户查看的流派：

top_rated_movies <- movies_df[as.numeric(names(inputUser)[order(inputUser,decreasing = TRUE)]),]$Title
top_rated_genres <- movies_df[as.numeric(names(inputUser)[order(inputUser,decreasing = TRUE)]),]$Genres
top_rated_genres <- as.data.frame(top_rated_genres,stringsAsFactors=F)
top_rated_genres$count <- 1
top_rated_genres <- aggregate(count~top_rated_genres,FUN=sum,data=top_rated_genres)
top_rated_genres <- top_rated_genres[with(top_rated_genres, order(-count)), ]
top_rated_genres$top_rated_genres <- factor(top_rated_genres$top_rated_genres, levels = top_rated_genres$top_rated_genres)
ggplot(top_rated_genres[top_rated_genres$count>1,],aes(x=top_rated_genres,y=count))+
geom_bar(stat="identity")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs(x="Genres",y="count",)+
theme(plot.title = element_text(hjust = 0.5))

重建输入向量以获取所有流派/电影的推荐百分比：

hh0 = tf$nn$sigmoid(tf$matmul(v0, W) + hb)
vv1 = tf$nn$sigmoid(tf$matmul(hh0, tf$transpose(W)) + vb)
feed = sess$run(hh0, feed_dict=dict( v0= inputUser, W= prv_w, hb= prv_hb))
rec = sess$run(vv1, feed_dict=dict( hh0= feed, W= prv_w, vb= prv_vb))
names(rec) <- movies_df$id_order

绘制顶级推荐流派：

top_recom_genres <- movies_df[as.numeric(names(rec)[order(rec,decreasing = TRUE)]),]$Genres
top_recom_genres <- as.data.frame(top_recom_genres,stringsAsFactors=F)
top_recom_genres$count <- 1
top_recom_genres <- aggregate(count~top_recom_genres,FUN=sum,data=top_recom_genres)
top_recom_genres <- top_recom_genres[with(top_recom_genres, order(-count)), ]
top_recom_genres$top_recom_genres <- factor(top_recom_genres$top_recom_genres, levels = top_recom_genres$top_recom_genres)
ggplot(top_recom_genres[top_recom_genres$count>20,],aes(x=top_recom_genres,y=count))+
geom_bar(stat="identity")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs(x="Genres",y="count",)+
theme(plot.title = element_text(hjust = 0.5))

查找前 10 个推荐电影：

top_recom_movies <- movies_df[as.numeric(names(rec)[order(rec,decreasing = TRUE)]),]$Title[1:10]

以下图片展示了前 10 个推荐电影：

设置深度置信网络（DBN）

深度置信网络（DBN）是一种深度神经网络（DNN），由多个隐藏层（或潜在变量）组成。在这里，连接只存在于各层之间，而层内的节点之间没有连接。DBN 可以作为无监督模型和有监督模型进行训练。

无监督模型用于去噪并重建输入，而有监督模型（经过预训练后）用于执行分类。由于每一层中的节点之间没有连接，DBN 可以看作是一个由无监督 RBM 或自编码器组成的集合，每个隐藏层作为其后续连接隐藏层的可见层。

这种堆叠 RBM 通过在所有层之间应用 CD，增强了输入重建的性能，从实际的输入训练层开始，直到最后的隐藏（或潜在）层。

DBN 是一种图模型，以贪婪的方式训练堆叠的 RBM。它们的网络倾向于使用输入特征向量i和隐藏层*h[1,2....m]*之间的联合分布来学习深层次的层次表示：

这里，i = h[0]；P(h[k-1]|h[k])是 RBM 第k层的隐藏层上重构可见单元的条件分布；*P(h[m-1],h[m])*是 DBN 最终 RBM 层的隐藏单元和可见单元（重构）的联合分布。下图展示了一个具有四个隐藏层的 DBN，其中W表示权重矩阵：

DBN 还可以增强 DNN 的鲁棒性。DNN 在实现反向传播时面临局部优化问题。在错误表面有许多低谷的情况下，反向传播导致梯度下降发生在局部深谷中（而不是全局深谷）。而 DBN 则对输入特征进行预训练，这有助于优化指向全局最深的低谷，然后使用反向传播执行梯度下降，以逐步最小化误差率。

训练三个 RBM 的堆叠：在本配方中，我们将使用三个堆叠的 RBM 训练一个 DBN，其中第一个隐藏层有 900 个节点，第二个隐藏层有 500 个节点，第三个隐藏层有 300 个节点。

准备好

本节提供了 TensorFlow 的要求。

数据集已加载并设置完毕
使用以下脚本加载TensorFlow包：

require(tensorflow)

怎么做...

本配方涵盖了设置深度置信网络（DBM）的步骤：

将每个隐藏层的节点数定义为一个向量：

RBM_hidden_sizes = c(900, 500 , 300 )

利用为伯努利分布输入设置限制玻尔兹曼机配方中展示的代码，生成一个 RBM 函数，使用以下输入和输出参数：

这是设置 RBM 的函数：

RBM <- function(input_data, num_input, num_output, epochs = 5, alpha = 0.1, batchsize=100){
# Placeholder variables
vb <- tf$placeholder(tf$float32, shape = shape(num_input))
hb <- tf$placeholder(tf$float32, shape = shape(num_output))
W <- tf$placeholder(tf$float32, shape = shape(num_input, num_output))
# Phase 1 : Forward Phase
X = tf$placeholder(tf$float32, shape=shape(NULL, num_input))
prob_h0= tf$nn$sigmoid(tf$matmul(X, W) + hb) #probabilities of the hidden units
h0 = tf$nn$relu(tf$sign(prob_h0 - tf$random_uniform(tf$shape(prob_h0)))) #sample_h_given_X
# Phase 2 : Backward Phase
prob_v1 = tf$nn$sigmoid(tf$matmul(h0, tf$transpose(W)) + vb)
v1 = tf$nn$relu(tf$sign(prob_v1 - tf$random_uniform(tf$shape(prob_v1))))
h1 = tf$nn$sigmoid(tf$matmul(v1, W) + hb)
# calculate gradients
w_pos_grad = tf$matmul(tf$transpose(X), h0)
w_neg_grad = tf$matmul(tf$transpose(v1), h1)
CD = (w_pos_grad - w_neg_grad) / tf$to_float(tf$shape(X)[0])
update_w = W + alpha * CD
update_vb = vb + alpha * tf$reduce_mean(X - v1)
update_hb = hb + alpha * tf$reduce_mean(h0 - h1)
# Objective function
err = tf$reduce_mean(tf$square(X - v1))
# Initialize variables
cur_w = tf$Variable(tf$zeros(shape = shape(num_input, num_output), dtype=tf$float32))
cur_vb = tf$Variable(tf$zeros(shape = shape(num_input), dtype=tf$float32))
cur_hb = tf$Variable(tf$zeros(shape = shape(num_output), dtype=tf$float32))
prv_w = tf$Variable(tf$random_normal(shape=shape(num_input, num_output), stddev=0.01, dtype=tf$float32))
prv_vb = tf$Variable(tf$zeros(shape = shape(num_input), dtype=tf$float32))
prv_hb = tf$Variable(tf$zeros(shape = shape(num_output), dtype=tf$float32))
# Start tensorflow session
sess$run(tf$global_variables_initializer())
output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(X=input_data,
W = prv_w$eval(),
vb = prv_vb$eval(),
hb = prv_hb$eval()))
prv_w <- output[[1]]
prv_vb <- output[[2]]
prv_hb <- output[[3]]
sess$run(err, feed_dict=dict(X= input_data, W= prv_w, vb= prv_vb, hb= prv_hb))
errors <- weights <- list()
u=1
for(ep in 1:epochs){
for(i in seq(0,(dim(input_data)[1]-batchsize),batchsize)){
batchX <- input_data[(i+1):(i+batchsize),]
output <- sess$run(list(update_w, update_vb, update_hb), feed_dict = dict(X=batchX,
W = prv_w,
vb = prv_vb,
hb = prv_hb))
prv_w <- output[[1]]
prv_vb <- output[[2]]
prv_hb <- output[[3]]
if(i%%10000 == 0){
errors[[u]] <- sess$run(err, feed_dict=dict(X= batchX, W= prv_w, vb= prv_vb, hb= prv_hb))
weights[[u]] <- output[[1]]
u=u+1
cat(i , " : ")
}
}
cat("epoch :", ep, " : reconstruction error : ", errors[length(errors)][[1]],"\n")
}
w <- prv_w
vb <- prv_vb
hb <- prv_hb
# Get the output
input_X = tf$constant(input_data)
ph_w = tf$constant(w)
ph_hb = tf$constant(hb)
out = tf$nn$sigmoid(tf$matmul(input_X, ph_w) + ph_hb)
sess$run(tf$global_variables_initializer())
return(list(output_data = sess$run(out),
error_list=errors,
weight_list=weights,
weight_final=w,
bias_final=hb))
}

按顺序训练三种不同类型的隐藏节点。换句话说，首先训练包含 900 个隐藏节点的 RBM1，然后将 RBM1 的输出作为输入，训练包含 500 个隐藏节点的 RBM2，再将 RBM2 的输出作为输入，训练包含 300 个隐藏节点的 RBM3。将所有三个 RBM 的输出存储为一个列表，RBM_output：

inpX = trainX
RBM_output <- list()
for(i in 1:length(RBM_hidden_sizes)){
size <- RBM_hidden_sizes[i]
# Train the RBM
RBM_output[[i]] <- RBM(input_data= inpX,
num_input= ncol(trainX),
num_output=size,
epochs = 5,
alpha = 0.1,
batchsize=100)
# Update the input data
inpX <- RBM_output[[i]]$output_data
# Update the input_size
num_input = size
cat("completed size :", size,"\n")
}

创建一个包含三个隐藏层批次误差的数据框：

error_df <- data.frame("error"=c(unlist(RBM_output[[1]]$error_list),unlist(RBM_output[[2]]$error_list),unlist(RBM_output[[3]]$error_list)),
"batches"=c(rep(seq(1:length(unlist(RBM_output[[1]]$error_list))),times=3)),
"hidden_layer"=c(rep(c(1,2,3),each=length(unlist(RBM_output[[1]]$error_list)))),
stringsAsFactors = FALSE)

绘制重构均方误差：

plot(error ~ batches,
xlab = "# of batches",
ylab = "Reconstruction Error",
pch = c(1, 7, 16)[hidden_layer],
main = "Stacked RBM-Reconstruction MSE plot",
data = error_df)
legend('topright',
c("H1_900","H2_500","H3_300"),
pch = c(1, 7, 16))

它是如何工作的...

评估训练三个堆叠 RBM 的性能：在这里，我们将为每个 RBM 运行五个时期（或迭代）。每个时期将执行批量优化（批大小 = 100）。在每个批次中，计算 CD 并相应地更新权重和偏置。

为了跟踪优化过程，在每批 10,000 行数据后计算均方误差。下图显示了为三种 RBM 分别计算的 30 个批次的均方重构误差下降趋势：

实现前馈反向传播神经网络

在这个配方中，我们将实现一个带有反向传播的柳树神经网络。神经网络的输入是第三个（或最后一个）RBM 的结果。换句话说，重构的原始数据（trainX）实际上被用来训练神经网络，作为一个监督分类器来识别（10）个数字。反向传播技术用于进一步调整分类性能。

准备工作

本节提供了 TensorFlow 的要求。

数据集被加载并设置完毕
TensorFlow包已经设置并加载

如何做…

本节涵盖了设置前馈反向传播神经网络的步骤：

让我们将神经网络的输入参数定义为函数参数。下表描述了每个参数：

神经网络函数将具有如下脚本中所示的结构：

NN_train <- function(Xdata,Ydata,Xtestdata,Ytestdata,input_size,
learning_rate=0.1,momentum = 0.1,epochs=10,
batchsize=100,rbm_list,dbn_sizes){
library(stringi)
*## insert all the codes mentioned in next 11 points* }

初始化一个长度为 4 的权重和偏置列表，第一个是一个 784 x 900 的标准差为 0.01 的正态分布张量，第二个是 900 x 500，第三个是 500 x 300，第四个是 300 x 10：

weight_list <- list()
bias_list <- list()
# Initialize variables
for(size in c(dbn_sizes,ncol(Ydata))){
#Initialize weights through a random uniform distribution
weight_list <- c(weight_list,tf$random_normal(shape=shape(input_size, size), stddev=0.01, dtype=tf$float32))
#Initialize bias as zeroes
bias_list <- c(bias_list, tf$zeros(shape = shape(size), dtype=tf$float32))
input_size = size
}

检查堆叠的 RBM 的结果是否符合dbn_sizes参数中提到的隐藏层尺寸：

#Check if expected dbn_sizes are correct
if(length(dbn_sizes)!=length(rbm_list)){
stop("number of hidden dbn_sizes not equal to number of rbm outputs generated")
# check if expected sized are correct
for(i in 1:length(dbn_sizes)){
if(dbn_sizes[i] != dbn_sizes[i])
stop("Number of hidden dbn_sizes do not match")
}
}

现在，将权重和偏置放置在weight_list和bias_list中的合适位置：

for(i in 1:length(dbn_sizes)){
weight_list[[i]] <- rbm_list[[i]]$weight_final
bias_list[[i]] <- rbm_list[[i]]$bias_final
}

创建输入和输出数据的占位符：

input <- tf$placeholder(tf$float32, shape = shape(NULL,ncol(Xdata)))
output <- tf$placeholder(tf$float32, shape = shape(NULL,ncol(Ydata)))

现在，使用从堆叠 RBM 中获得的权重和偏置来重构输入数据，并将每个 RBM 的重构数据存储在列表input_sub中：

input_sub <- list()
weight <- list()
bias <- list()
for(i in 1:(length(dbn_sizes)+1)){
weight[[i]] <- tf$cast(tf$Variable(weight_list[[i]]),tf$float32)
bias[[i]] <- tf$cast(tf$Variable(bias_list[[i]]),tf$float32)
}
input_sub[[1]] <- tf$nn$sigmoid(tf$matmul(input, weight[[1]]) + bias[[1]])
for(i in 2:(length(dbn_sizes)+1)){
input_sub[[i]] <- tf$nn$sigmoid(tf$matmul(input_sub[[i-1]], weight[[i]]) + bias[[i]])
}

定义成本函数——即预测与实际数字之间的均方误差：

cost = tf$reduce_mean(tf$square(input_sub[[length(input_sub)]] - output))

实现反向传播以最小化成本：

train_op <- tf$train$MomentumOptimizer(learning_rate, momentum)$minimize(cost)

生成预测结果：

predict_op = tf$argmax(input_sub[[length(input_sub)]],axis=tf$cast(1.0,tf$int32))

执行训练迭代：

train_accuracy <- c()
test_accuracy <- c()
for(ep in 1:epochs){
for(i in seq(0,(dim(Xdata)[1]-batchsize),batchsize)){
batchX <- Xdata[(i+1):(i+batchsize),]
batchY <- Ydata[(i+1):(i+batchsize),]
#Run the training operation on the input data
sess$run(train_op,feed_dict=dict(input = batchX,
output = batchY))
}
for(j in 1:(length(dbn_sizes)+1)){
# Retrieve weights and biases
weight_list[[j]] <- sess$run(weight[[j]])
bias_list[[j]] <- sess$ run(bias[[j]])
}
train_result <- sess$run(predict_op, feed_dict = dict(input=Xdata, output=Ydata))+1
train_actual <- as.numeric(stringi::stri_sub(colnames(as.data.frame(Ydata))[max.col(as.data.frame(Ydata),ties.method="first")],2))
test_result <- sess$run(predict_op, feed_dict = dict(input=Xtestdata, output=Ytestdata))+1
test_actual <- as.numeric(stringi::stri_sub(colnames(as.data.frame(Ytestdata))[max.col(as.data.frame(Ytestdata),ties.method="first")],2))
train_accuracy <- c(train_accuracy,mean(train_actual==train_result))
test_accuracy <- c(test_accuracy,mean(test_actual==test_result))
cat("epoch:", ep, " Train Accuracy: ",train_accuracy[ep]," Test Accuracy : ",test_accuracy[ep],"\n")
}

最后，返回四个结果的列表，分别是训练准确度（train_accuracy）、测试准确度（test_accuracy）、每次迭代生成的权重矩阵列表（weight_list）和每次迭代生成的偏置向量列表（bias_list）：

return(list(train_accuracy=train_accuracy,
test_accuracy=test_accuracy,
weight_list=weight_list,
bias_list=bias_list))

对定义的神经网络进行训练迭代：

NN_results <- NN_train(Xdata=trainX,
Ydata=trainY,
Xtestdata=testX,
Ytestdata=testY,
input_size=ncol(trainX),
rbm_list=RBM_output,
dbn_sizes = RBM_hidden_sizes)

以下代码用于绘制训练和测试准确度：

accuracy_df <- data.frame("accuracy"=c(NN_results$train_accuracy,NN_results$test_accuracy),
"epochs"=c(rep(1:10,times=2)),
"datatype"=c(rep(c(1,2),each=10)),
stringsAsFactors = FALSE)
plot(accuracy ~ epochs,
xlab = "# of epochs",
ylab = "Accuracy in %",
pch = c(16, 1)[datatype],
main = "Neural Network - Accuracy in %",
data = accuracy_df)
legend('bottomright',
c("train","test"),
pch = c( 16, 1))

它是如何工作的…

评估神经网络的训练和测试性能：下图显示了在训练神经网络时观察到的训练和测试准确度的增长趋势：

设置深度限制玻尔兹曼机

与 DBNs 不同，深度限制玻尔兹曼机 (DRBM) 是由互联的隐藏层组成的无向网络，能够学习这些连接上的联合概率。在当前的设置中，每次迭代后，能见度和隐藏变量会从偏置偏移向量中减去，进行中心化处理。研究表明，中心化优化了 DRBM 的性能，并且与传统的 RBM 相比，可以达到更高的对数似然值。

准备就绪

本节提供了设置 DRBM 的要求：

已加载并设置 MNIST 数据集：
已设置并加载 tensorflow 包：

如何操作...

本节详细介绍了如何在 R 中使用 TensorFlow 设置 DRBM 模型的步骤：

定义 DRBM 的参数：

learning_rate = 0.005
momentum = 0.005
minbatch_size = 25
hidden_layers = c(400,100)
biases = list(-1,-1)

使用双曲正切定义一个 sigmoid 函数 [(log(1+x) -log(1-x))/2]：

arcsigm <- function(x){
return(atanh((2*x)-1)*2)
}

使用双曲正切定义一个 sigmoid 函数 [(e^x-e^(-x))/(e^x+e^(-x))]：

sigm <- function(x){
return(tanh((x/2)+1)/2)
}

定义一个 binarize 函数，返回一个二值矩阵（0,1）：

binarize <- function(x){
# truncated rnorm
trnrom <- function(n, mean, sd, minval = -Inf, maxval = Inf){
qnorm(runif(n, pnorm(minval, mean, sd), pnorm(maxval, mean, sd)), mean, sd)
}
return((x > matrix( trnrom(n=nrow(x)*ncol(x),mean=0,sd=1,minval=0,maxval=1), nrow(x), ncol(x)))*1)
}

定义一个 re_construct 函数，返回一个像素矩阵：

re_construct <- function(x){
x = x - min(x) + 1e-9
x = x / (max(x) + 1e-9)
return(x*255)
}

定义一个函数来执行给定层的 gibbs 激活：

gibbs <- function(X,l,initials){
if(l>1){
bu <- (X[l-1][[1]] - matrix(rep(initials$param_O[[l-1]],minbatch_size),minbatch_size,byrow=TRUE)) %*%
initials$param_W[l-1][[1]]
} else {
bu <- 0
}
if((l+1) < length(X)){
td <- (X[l+1][[1]] - matrix(rep(initials$param_O[[l+1]],minbatch_size),minbatch_size,byrow=TRUE))%*%
t(initials$param_W[l][[1]])
} else {
td <- 0
}
X[[l]] <- binarize(sigm(bu+td+matrix(rep(initials$param_B[[l]],minbatch_size),minbatch_size,byrow=TRUE)))
return(X[[l]])
}

定义一个函数来执行偏置向量的重参数化：

reparamBias <- function(X,l,initials){
if(l>1){
bu <- colMeans((X[[l-1]] - matrix(rep(initials$param_O[[l-1]],minbatch_size),minbatch_size,byrow=TRUE))%*%
initials$param_W[[l-1]])
} else {
bu <- 0
}
if((l+1) < length(X)){
td <- colMeans((X[[l+1]] - matrix(rep(initials$param_O[[l+1]],minbatch_size),minbatch_size,byrow=TRUE))%*%
t(initials$param_W[[l]]))
} else {
td <- 0
}
initials$param_B[[l]] <- (1-momentum)*initials$param_B[[l]] + momentum*(initials$param_B[[l]] + bu + td)
return(initials$param_B[[l]])
}

定义一个函数来执行偏置偏移向量的重参数化：

reparamO <- function(X,l,initials){
initials$param_O[[l]] <- colMeans((1-momentum)*matrix(rep(initials$param_O[[l]],minbatch_size),minbatch_size,byrow=TRUE) + momentum*(X[[l]]))
return(initials$param_O[[l]])
}

定义一个函数来初始化权重、偏置、偏置偏移和输入矩阵：

DRBM_initialize <- function(layers,bias_list){
# Initialize model parameters and particles
param_W <- list()
for(i in 1:(length(layers)-1)){
param_W[[i]] <- matrix(0L, nrow=layers[i], ncol=layers[i+1])
}
param_B <- list()
for(i in 1:length(layers)){
param_B[[i]] <- matrix(0L, nrow=layers[i], ncol=1) + bias_list[[i]]
}
param_O <- list()
for(i in 1:length(param_B)){
param_O[[i]] <- sigm(param_B[[i]])
}
param_X <- list()
for(i in 1:length(layers)){
param_X[[i]] <- matrix(0L, nrow=minbatch_size, ncol=layers[i]) + matrix(rep(param_O[[i]],minbatch_size),minbatch_size,byrow=TRUE)
}
return(list(param_W=param_W,param_B=param_B,param_O=param_O,param_X=param_X))
}

使用前面配方中引入的 MNIST 训练数据（trainX）。通过将 trainX 数据除以 255，进行标准化处理：

X <- trainX/255

生成初始的权重矩阵、偏置向量、偏置偏移向量和输入矩阵：

layers <- c(784,hidden_layers)
bias_list <- list(arcsigm(pmax(colMeans(X),0.001)),biases[[1]],biases[[2]])
initials <-DRBM_initialize(layers,bias_list)

对输入数据 X 进行子集化（minbatch_size）：

batchX <- X[sample(nrow(X))[1:minbatch_size],]

执行 1,000 次迭代。在每次迭代中，更新初始权重和偏置 100 次，并绘制权重矩阵的图像：

for(iter in 1:1000){

# Perform some learnings
for(j in 1:100){
# Initialize a data particle
dat <- list()
dat[[1]] <- binarize(batchX)
for(l in 2:length(initials$param_X)){
dat[[l]] <- initials$param_X[l][[1]]*0 + matrix(rep(initials$param_O[l][[1]],minbatch_size),minbatch_size,byrow=TRUE)
}

# Alternate gibbs sampler on data and free particles
for(l in rep(c(seq(2,length(initials$param_X),2), seq(3,length(initials$param_X),2)),5)){
dat[[l]] <- gibbs(dat,l,initials)
}
for(l in rep(c(seq(2,length(initials$param_X),2), seq(1,length(initials$param_X),2)),1)){
initials$param_X[[l]] <- gibbs(initials$param_X,l,initials)
}

# Parameter update
for(i in 1:length(initials$param_W)){
initials$param_W[[i]] <- initials$param_W[[i]] + (learning_rate*((t(dat[[i]] - matrix(rep(initials$param_O[i][[1]],minbatch_size),minbatch_size,byrow=TRUE)) %*%
(dat[[i+1]] - matrix(rep(initials$param_O[i+1][[1]],minbatch_size),minbatch_size,byrow=TRUE))) -
(t(initials$param_X[[i]] - matrix(rep(initials$param_O[i][[1]],minbatch_size),minbatch_size,byrow=TRUE)) %*%
(initials$param_X[[i+1]] - matrix(rep(initials$param_O[i+1][[1]],minbatch_size),minbatch_size,byrow=TRUE))))/nrow(batchX))
}

for(i in 1:length(initials$param_B)){
initials$param_B[[i]] <- colMeans(matrix(rep(initials$param_B[[i]],minbatch_size),minbatch_size,byrow=TRUE) + (learning_rate*(dat[[i]] - initials$param_X[[i]])))
}

# Reparameterization
for(l in 1:length(initials$param_B)){
initials$param_B[[l]] <- reparamBias(dat,l,initials)
}
for(l in 1:length(initials$param_O)){
initials$param_O[[l]] <- reparamO(dat,l,initials)
}
}

# Generate necessary outputs
cat("Iteration:",iter," ","Mean of W of VL-HL1:",mean(initials$param_W[[1]])," ","Mean of W of HL1-HL2:",mean(initials$param_W[[2]]) ,"\n")
cat("Iteration:",iter," ","SDev of W of VL-HL1:",sd(initials$param_W[[1]])," ","SDev of W of HL1-HL2:",sd(initials$param_W[[2]]) ,"\n")

# Plot weight matrices
W=diag(nrow(initials$param_W[[1]]))
for(l in 1:length(initials$param_W)){
W = W %*% initials$param_W[[l]]
m = dim(W)[2] * 0.05
w1_arr <- matrix(0,28*m,28*m)
i=1
for(k in 1:m){
for(j in 1:28){
vec <- c(W[(28*j-28+1):(28*j),(k*m-m+1):(k*m)])
w1_arr[i,] <- vec
i=i+1
}
}
w1_arr = re_construct(w1_arr)
w1_arr <- floor(w1_arr)
image(w1_arr,axes = TRUE, col = grey(seq(0, 1, length = 256)))
}
}

它是如何工作的...

由于前面的 DRBM 使用了两个隐藏层进行训练，我们生成了两个权重矩阵。第一个权重矩阵定义了可见层和第一个隐藏层之间的连接。第二个权重矩阵定义了第一个隐藏层和第二个隐藏层之间的连接。下图显示了第一个权重矩阵的像素图像：

下图显示了第二个权重矩阵的第二个像素图像：