基于无监督学习的自编码器实现

无监督学习介绍及部分函数实现

1. 内容简介

   目前许多有监督学习算法，如SVM，DNN或是boosting，决策树等，都在工业界分类或决策任务上取得了不错的成果。但是这些有监督学习需要大量带标签的数据，对数据进行上标签又是一个需耗费人力与时间的任务。有许多数据都是不带标签的，因此我们可利用无监督学习对其进行聚类或特征提取。利用无监督学习得到的特征结果也可应用到带标记数据较少的有监督学习任务中，提高其分类性能。

   这门课程，我们将介绍基于无监督学习的一个简单应用——

   自编码器 ，一个可用于降维或者提取特征的神经网络。

2. 实验原理

   无监督学习的目的是利用无标记数据推断出数据内部隐藏的结构特征。与有监督学习从有标记的数据中学到一个有正确答案的模型不同，由于无监督学习的数据是未经标记的，因此在学习过程中就没有所谓误差或者指导信号去评估一个可能的解决方案。

我们可以从下图中看出无监督学习与有监督学习的区别：

与有监督学习有明确的目标输出不同，无监督学习没有给定目标输出，而是去提取数据本身的静态结构特征，如在上图中，有监督学习的目的是可对输入的图形进行分类，区分这些图形是什么形状的，但是无监督学习开始训练前，对于机器来说这些图形没有什么区别，但在经过训练后，我们能够获得这些数据的结构特征，根据这些图形各自的特点，根据图形的相似点对其进行聚类。要注意的是，在无监督学习之后，我们也不知道这些图形分别属于什么几何形状，仅是根据物以类聚的概念，将其划分为几个不同的集群。

如上描述的，就是一个常见的无监督学习 ——数据聚类。在人工神经网络中，自我组织映射和适应性共振理论 也是常用的无监督学习，也有许多网络是利用无监督学习对数据进行特征提取，如在图像识别的应用中，可用无监督学习提取图片中物体的边缘特征。

自编码器

在该门课中我们会介绍，如何利用目标输出与网络实际输出的误差进行反向传播，调节每一层间的权值矩阵。自编码器是基于无监督学习为对数据进行更加有效的编码的应用， 自编码器将目标输出等同于输入，同样也是使用误差反向传播机制进行权值调整。自编码器 是一个数据降维压缩算法，算法训练的结果是拟合一个恒等函数

$$h_{w,b}(x)=x$$ , 即通过自编码器的训练后， $$decoder(encoder(x))=x$$
目前自编码器主要有 数据降噪和为 数据可视化而降维两个方面的应用。这里我们设置的自编码器是具有一层隐藏节点的三层全连接神经网络，实现的是为数据可视化而降维这方面的应用，上图中展示的就是对数字图像 ‘3’ 降维后载重现的效果。

本次课程中，我们实现的是基于minist手写体数字数据库中的500张图片，每个数字 (0-9)都有其对应的50张手写体图片，图片的维度为28 * 28，在我们设置的网络中，手写体数字图片是作为网络的外部输入结点，内部输入节点设置为1个，隐藏层结点为32个，输出结点与外部输入节点个数一致。 即使通过我们的训练的自编码器，可将原始数据的维度从 784784 降维至 3232。

在我们的自编码器的三层神经网络中，根据网络的输入前向计算网络中所有节点的值时，所用的激活函数为sigmoid函数：

故在前向计算过程中，神经网络中隐藏层节点的值为:

$$h_i=f(z^h_i)$$
$$z^h_i=\sum_{j=1}^nx_jw_{j,i}+aw_{n+1,j}$$
输出层的节点值为：
$$\hat{x_i}=f(\hat{z_i})$$
$$\hat{z_i}=\sum_{j=1}^mw_{j,i}h_j$$
自编码器网络计算实际输出与输入值之间的误差函数（Cost function）CC为：
$$C=1/2*\sum_{i=1}^n(\hat{x_i}-y)^2$$
故在误差反向传播时，第ll层的δ值为：

这里所实现的自编码器并没有对编码层（隐藏层）作稀疏限制，若加上稀疏限制使编码层具有稀疏性，则在误差代价函数和误差反向传播就要做相应的修改。

稀疏性：只有很少的几个非零元素或只有很少的几个远大于零的元素

---

自编码器实现

1. 实验过程

   数据预处理
   在实验过程中，我们用一个列向量表示一个图片样本，如数字图片维度为28*2828∗28时，我们用784*1784∗1的列向量表示其数据，这个过程称为向量化数据。由于输入灰度图片像素值范围在0~255，为加快网络的训练收敛，且原始图像数据由于是0-255的unit数据为便于python里的浮点型运算，我们将图像数据归一化至0-1区间。
   这里我们采用的归一化方法为最小-最大标准化（Min-Max Nomalization），转换函数如下：
   

   $$x^*=(x-min)/(max-min)$$

   由于这里图像数据的最大值为255，最小值为0，故在实验过程中，我们归一化时直接将原始数据除以255即可。
   

自编码器训练过程
众所周知，机器学习的三大类可分为：

1. 监督学习
2. 无监督学习
3. 强化学习

其实还有一种属于这两种学习的折中，称为半监督学习，同时利用无标记和有标记的数据进行训练。

若按照学习过程(权值的更新方法)中，一次使用多少个样本对权值进行更新划分，可划分成：

1. 离线学习 / 批梯度下降
2. 在线学习 / 随机梯度下降
3. 最小批梯度下降

批梯度下降是在所有的数据样本训练完毕之后，累加所有样本对权值的更新值，对权值进行一次更新，在训练过程中当有新的样本出现时，我们的权值也没有进行及时的更新，故也称为离线学习。 随机梯度下降是每次只用一个样本对权值进行更新，所以当新的样本出现时，我们也能及时更新权值，故也称为在线学习。

批梯度下降算法需要所有的样本进行一次权值更新，若是数据集中有很多样本时，每次更新说花费的时间就很多。而随机梯度下降虽比批梯度下降算法更行权值速度快，但是由于更新得过于频繁，且单个样本所得到的代价函数比用所有样本得到的代价函数更具有随机性，可能无法找到逼近于我们想要的恒等函数$h_{w,b}(x)=x$的最优解。故我们采取了一种折中的学习过程，mini-batch learning最小批训练。即每次训练时，只用所有数据中的一部分数据对权值进行更新。

在我们的训练过程中，数据样本集大小为500，每次更新权值时使用100张图像对其进行更新。这里我们设置的总训练次数（即迭代次数）为500次，也可根据代价函数 $C =\frac{1}{2}\sum_{i=1}^{n}(\hat{x}_i - x_i)^2$
​​ ，当根据代价函数计算所得的误差小于某个 任意小的$\epsilon, \epsilon>0ϵ,ϵ>0$时停止训练。

代码实现

#-*- coding=utf8 -*-

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
import random


def main():

    trainData = scio.loadmat(r'D:\python_learning_dataes\trainData.mat')

    unlabeled_data = trainData['trainData']
    unlabeled_data = unlabeled_data[:,:] / 255.


    # define the learning parameters
    alpha = 5 # learning rate
    max_epoch = 500 # the learning epoch
    mini_batch = 100 # used for the batch learning
    imgSize = 784 # image size of a digit image in the minist dataset

    # the network structure
    layer_struc = [[imgSize, 1],
                   [0, 32],
                   [0, imgSize]]
    layer_num = 3

    # initialize weights
    w = []
    for l in range(layer_num-1):
        w.append(np.random.randn(layer_struc[l+1][1],sum(layer_struc[l])))

    dataset_size = 500  # the number of the images in the training data set

    # define the internal input of the network
    X = []
    X.append(np.array(unlabeled_data[:,:]))
    X.append(np.zeros((0,dataset_size)))
    X.append(np.zeros((0,dataset_size)))

    delta = []
    for l in range(layer_num):
        delta.append([])

    # define the display parameters 
    nRow = max_epoch / 100 + 1
    nColumn = 10   # display 10 images in each row 
    eachDigitNum = 50  # 50 instancese corresponding to each digit in the training set

    # display the original digit in the first row

    for iImg in range(nColumn):
        ax = plt.subplot(nRow, nColumn, iImg+1)
        plt.imshow(unlabeled_data[:,eachDigitNum * iImg + 1].reshape((28,28)).T, cmap= plt.cm.gray)

        if iImg == 0:
            plt.ylabel('Original Images',rotation=90)

        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

    # unsupervised training
    count = 0 # count the iteration
    print('Autoencoder training start..')
    for iter in range(max_epoch):

        # define the shuffle index
        ind = list(range(dataset_size))
        random.shuffle(ind)

        a = []
        z = []
        z.append([])
        for i in range(int(np.ceil(dataset_size / mini_batch))):
            a.append(np.zeros((layer_struc[0][1], mini_batch)))
            x = []
            for l in range(layer_num):
                x.append( X[l][:,ind[i*mini_batch : min((i+1)*mini_batch, dataset_size)]])

            y = unlabeled_data[:,ind[i*mini_batch:min((i+1)*mini_batch,dataset_size)]]
            for l in range(layer_num-1):
                a.append([])
                z.append([])
                a[l+1],z[l+1] = feedforward(w[l],a[l],x[l])


            delta[layer_num-1] = np.array(a[layer_num-1] - y) * np.array(a[layer_num-1])
            delta[layer_num-1] = delta[layer_num-1] * np.array(1-a[layer_num-1])

            for l in range(layer_num-2, 0, -1):
                delta[l] = backprop(w[l],z[l],delta[l+1])

            for l in range(layer_num-1):
                dw = np.dot(delta[l+1], np.concatenate((a[l],x[l]),axis=0).T) / mini_batch
                w[l] = w[l] - alpha * dw

        count = count + 1  



        # display reconstruction result 
        if np.mod(iter+1,100) == 0 :
            b = []
            b.append(np.zeros((layer_struc[0][1],dataset_size)))

            for l in range(layer_num-1):
                tempA, tempZ = feedforward(w[l], b[l], X[l])                
                b.append(tempA)

            for iImg in range(nColumn):
                ax = plt.subplot(nRow,nColumn, iImg + nColumn * (iter+1)/100 + 1)
                dis_result = b[layer_num-1][:,eachDigitNum * iImg + 1].reshape(28,28).T
                plt.imshow(dis_result,cmap= plt.cm.gray) 
                ax.get_xaxis().set_visible(False)
                ax.get_yaxis().set_visible(False)

            print('Learning epoch:', count, '/', max_epoch)

    plt.show()


# feedforward computing
def feedforward(w,a,x):
    f = lambda s: 1 / (1 + np.exp(-s)) 

    # concatenate the matrix a and x , and multiply with weight matrix
    w = np.array(w)
    temp = np.array(np.concatenate((a,x),axis=0))
    z_next = np.dot(w , temp)

    return f(z_next), z_next

# backpropagation
def backprop(w,z,delta_next):

    # sigmoid function
    f = lambda s: np.array(1 / (1 + np.exp(-s)))

    # the Derivative of sigmoid function
    df = lambda s: f(s) * (1 - f(s))


    delta = df(z) * np.dot(w.T,delta_next)    

    return delta

if __name__ == '__main__':
    main()