theano-多分类逻辑回归代码解析

Theano的多分类逻辑回归代码解析

这份代码主要借助theano框架用多分类的逻辑回归模型做手写体数字识别的问题，优化模型时采用的是随机
梯度下降的算法。

下面是多分类逻辑回归的模型公式，它的参数是权值矩阵W和偏值向量b,最后公式(1)中的P得出的是属
于这一类的概率，而整个模型的输出是所有$P(Y=i|x, W,b)$中最大的一个，例如一个10分类问题，最终计
算出10个P，最大的那个P所对应的类别i(注意这里有多少类，i就有几个)就是最终的分类输出类别。

$$P(Y=i|x, W,b)=softmax_i(Wx+b)=\frac{e^{W_ix+b_i}}{\sum_je^{W_jx+b_j}}\qquad(1)$$
$$y_{pred}=argmax_iP(Y=i|x,W,b)\qquad(2)$$

1. 这部分是导入需要的模块

from __future__ import print_function

import six.moves.cPickle as pickle
import gzip
import os
import sys
import timeit
import numpy
import theano
import theano.tensor as T

__docformat__ = 'restructedtext en'

2. 这部分是Logistic回归的类定义

class LogisticRegression(object):

    def __init__(self, input, n_in, n_out):
        """ 初始化所有模型参数
        :type input: theano.tensor.TensorType
        :param input: symbolic variable that describes the input of the
                      architecture (one minibatch)

        :type n_in: int
        :param n_in: number of input units, the dimension of the space in
                     which the datapoints lie

        :type n_out: int
        :param n_out: number of output units, the dimension of the space in
                      which the labels lie
        """
        # 初始化权值矩阵为一个零矩阵
        self.W = theano.shared(
            value=numpy.zeros(
                (n_in, n_out),
                dtype=theano.config.floatX
            ),
            name='W',
            borrow=True
        )
        # 初始化偏值向量为零向量，大小为输出单元的大小，即类别大小
        self.b = theano.shared(
            value=numpy.zeros(
                (n_out,),
                dtype=theano.config.floatX
            ),
            name='b',
            borrow=True
        )

        # 根据模型输入计算输出概率，这里的p_y_given_x是一个向量，有多少个类别，这个向量就有多少个概率。
        self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)

        # 计算p_y_given_x向量中的最大值
        self.y_pred = T.argmax(self.p_y_given_x, axis=1)

        # 模型参数
        self.params = [self.W, self.b]

        # 模型输入
        self.input = input

现在初始化了模型输入，给定数据，就可以计算模型输出了，但是怎么计算这个数据分布下的模型的损失函数
那，我们采用最大似然估计的方法，下面的公式(3)就是似然函数，公式(4)就是损失函数，即损失函数其实
就是负对数似然函数。

$$\mathcal{L}{(\theta=\{W,b\},\mathcal{D})=\sum_{i=0}^{|\mathcal{D}|}}log(P(Y=y^ {(i)}|x^{(i)},W,b))\qquad(3)$$

$$\ell (\theta=\{W,b\},\mathcal{D})=-\mathcal{L}(\theta=\{W,b\},\mathcal{D})\qquad(4)$$

这里负对数似然函数的计算的代码比较难理解，虽然只有简单的一行代码，因为这里进行的都是符号计算，而
且theano是出了名的难调试，不好输出变量观察，所以这里仔细解释一个，为什么下面的那一行代码就可以
计算出对于给定minibatch数据的负对数似然函数的期望。这里y是训练数据集中标签数据的一个minibatch，
由于一开始使用的是loaddata函数导入的数据集，为了使用GPU计算，把所有数据转换成了TensorVariable
类型。

注：为讲解所需示例代码段，与整个程序无关
datasets = logistic_sgd.load_data('mnist.pkl.gz')
train_set_x,train_set_y = datasets[0]
y = train_set_y[0:10]
print(y.eval())
------output [0 4 1 9 2 1 3 1 4 3]

num = y.shape[0]
print(num.eval())
------output 10

a = T.arange(10)
print(a.eval())
------output [0 1 2 3 4 5 6 7 8 9]

从上面的代码段可以看出最终输出的是一个行向量，一共10个数据，表示这个minibatch取了10个样本，这
是10个样本的标签，每个标签是一个数字，代表的是类别，所以我们就知道了y.shape[0])其实代表的是样
本的个数，所以T.arange(y.shape[0])是一个大小为样本个数的有序向量。而p_y_given_x是似然函数，
每一个样本x经过p_y_given_x计算得到一个行向量，存放的是每一类的概率，而一个minibatch的x则输出
的p_y_given_x是一个矩阵，矩阵每行代表一个样本的输出概率，每列代表所有样本这一类的输出概率，这
里对p_y_given_x取了对数，则T.log(self.p_y_given_x)其实还是一个矩阵，我们假设用LP表示，那么
对于一个有10个样本的minibatch计算对数似然函数，下面那一行代码的结果就可以表示成
LP[[0 1 2 3 4 5 6 7 8 9],[0 4 1 9 2 1 3 1 4 3]]，注意这里LP的维度，行数对应样本个数，
列数对应类别个数，所以最后这个LP[]得到的正好是每个样本的输出概率中取到它对应的标签的那个概率，
最后得到的一个向量，这个向量存放的是所有样本最后计算出来的概率，然后把这些概率取均值，然后在取
负数，即得到负对数似然函数对于这个minibatch的期望，我们最终的目的是把这个期望最小化，也就是最小
化损失函数。

    def negative_log_likelihood(self, y):
        """ 计算负对数似然函数
        :type y: theano.tensor.TensorType
        :param y: corresponds to a vector that gives for each example the
                  correct label
        """
        return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])

下面的代码段是计算给定minibatch数据集上面，模型预测的错误率的。

    def errors(self, y):
        """计算最终模型的错误率
        :type y: theano.tensor.TensorType
        :param y: corresponds to a vector that gives for each example the
                  correct label
        """

        # 检查真实样本标签y和模型输出的y维度是否一致
        if y.ndim != self.y_pred.ndim:
            raise TypeError(
                'y should have the same shape as self.y_pred',
                ('y', y.type, 'y_pred', self.y_pred.type)
            )
        # 检查y的类型是否正确
        if y.dtype.startswith('int'):
            # T.neq操作输出的是一个包含0和1的向量，1表示预测错误，最后给这个向量取均值
            # 则得到了模型在这个minibatch上面的预测错误率
            return T.mean(T.neq(self.y_pred, y))
        else:
            raise NotImplementedError()

这段代码不做过多解释，因为和算法本身关系不是很大，是将数据集载入，并且转换成方便使用的类型

def load_data(dataset):
    ''' Loads the dataset

    :type dataset: string
    :param dataset: the path to the dataset (here MNIST)
    '''

    #############
    # LOAD DATA #
    #############

    # Download the MNIST dataset if it is not present
    data_dir, data_file = os.path.split(dataset)
    if data_dir == "" and not os.path.isfile(dataset):
        # Check if dataset is in the data directory.
        new_path = os.path.join(
            os.path.split(__file__)[0],
            "..",
            "data",
            dataset
        )
        if os.path.isfile(new_path) or data_file == 'mnist.pkl.gz':
            dataset = new_path

    if (not os.path.isfile(dataset)) and data_file == 'mnist.pkl.gz':
        from six.moves import urllib
        origin = (
            'http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz'
        )
        print('Downloading data from %s' % origin)
        urllib.request.urlretrieve(origin, dataset)

    print('... loading data')

    # Load the dataset
    with gzip.open(dataset, 'rb') as f:
        try:
            train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
        except:
            train_set, valid_set, test_set = pickle.load(f)
    # train_set, valid_set, test_set format: tuple(input, target)
    # input is a numpy.ndarray of 2 dimensions (a matrix)
    # where each row corresponds to an example. target is a
    # numpy.ndarray of 1 dimension (vector) that has the same length as
    # the number of rows in the input. It should give the target
    # to the example with the same index in the input.

    def shared_dataset(data_xy, borrow=True):
        """ Function that loads the dataset into shared variables

        The reason we store our dataset in shared variables is to allow
        Theano to copy it into the GPU memory (when code is run on GPU).
        Since copying data into the GPU is slow, copying a minibatch everytime
        is needed (the default behaviour if the data is not in a shared
        variable) would lead to a large decrease in performance.
        """
        data_x, data_y = data_xy
        shared_x = theano.shared(numpy.asarray(data_x,
                                               dtype=theano.config.floatX),
                                 borrow=borrow)
        shared_y = theano.shared(numpy.asarray(data_y,
                                               dtype=theano.config.floatX),
                                 borrow=borrow)
        # When storing data on the GPU it has to be stored as floats
        # therefore we will store the labels as ``floatX`` as well
        # (``shared_y`` does exactly that). But during our computations
        # we need them as ints (we use labels as index, and if they are
        # floats it doesn't make sense) therefore instead of returning
        # ``shared_y`` we will have to cast it to int. This little hack
        # lets ous get around this issue
        return shared_x, T.cast(shared_y, 'int32')

    test_set_x, test_set_y = shared_dataset(test_set)
    valid_set_x, valid_set_y = shared_dataset(valid_set)
    train_set_x, train_set_y = shared_dataset(train_set)

    rval = [(train_set_x, train_set_y), (valid_set_x, valid_set_y),
            (test_set_x, test_set_y)]
    return rval

这段代码也就是sgd_optimization_mnist函数是建立整个模型，并采用sgd进行优化的所有代码，下面代
码先载入数据集，然后从数据集中分割出来训练数据，测试数据和交叉验证数据。然后分别计算出针对给定
batch_size训练是需要迭代的次数，即整个数据集的个数除以batch_size，迭代完成一次就叫做一个epoch，
一般为了提高精度，会为epoch设置一个较大的数，即对整个数据集进行多次循环迭代。

def sgd_optimization_mnist(learning_rate=0.13, n_epochs=1000,
                           dataset='mnist.pkl.gz',
                           batch_size=600):
    ""
    :type learning_rate: float
    :param learning_rate: learning rate used (factor for the stochastic
                          gradient)

    :type n_epochs: int
    :param n_epochs: maximal number of epochs to run the optimizer

    :type dataset: string
    :param dataset: the path of the MNIST dataset file from
                 http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz

    """
    datasets = load_data(dataset)

    train_set_x, train_set_y = datasets[0]
    valid_set_x, valid_set_y = datasets[1]
    test_set_x, test_set_y = datasets[2]

    # compute number of minibatches for training, validation and testing
    n_train_batches = train_set_x.get_value(borrow=True).shape[0] // batch_size
    n_valid_batches = valid_set_x.get_value(borrow=True).shape[0] // batch_size
    n_test_batches = test_set_x.get_value(borrow=True).shape[0] // batch_size

下面这部分代码主要用于建立模型，主要定义几个内容，一个是计算模型正确率的函数，一个是计算验证数据
集正确率的函数，很重要的通过求损失函数的偏导数生成一个updates参数更新原则，然后再根据updates
建立一个train_model函数。

    ######################
    # BUILD ACTUAL MODEL #
    ######################
    print('... building the model')

    # allocate symbolic variables for the data
    index = T.lscalar()  # index to a [mini]batch

    # 生成符号变量 (x和y都表示minibatch)
    x = T.matrix('x')  # data, presented as rasterized images
    y = T.ivector('y')  # labels, presented as 1D vector of [int] labels

    # 每一个手写体数字大小是28乘以28,最终输出类别是10个
    classifier = LogisticRegression(input=x, n_in=28 * 28, n_out=10)

    # 损失函数是给定标签时的负对数似然函数
    cost = classifier.negative_log_likelihood(y)

    # 生成一个函数，givens是给定的一个minibatch大小的数据集，输入变量inputs表示选择的minibatch
    # 的索引，输出计算的是模型的错误率。
    test_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={
            x: test_set_x[index * batch_size: (index + 1) * batch_size],
            y: test_set_y[index * batch_size: (index + 1) * batch_size]
        }
    )

    # 这是一个交叉验证的函数，其实就是对于验证数据集，计算分类器模型的错误率的
    validate_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={
            x: valid_set_x[index * batch_size: (index + 1) * batch_size],
            y: valid_set_y[index * batch_size: (index + 1) * batch_size]
        }
    )

    # 计算损失函数cost关于模型参数theta = (W,b)的偏导数
    g_W = T.grad(cost=cost, wrt=classifier.W)
    g_b = T.grad(cost=cost, wrt=classifier.b)


    # 这里定义所有模型参数的更新方式，最后存放在一个列表list中，这个列表叫updates，列表中每个
    # 元素的类型是元组，每个元组里面是两个量，第一个量是参数，第二个量是这个参数的更新方式，即
    # (variable, update expression)这种形式。
    # specify how to update the parameters of the model as a list of
    # (variable, update expression) pairs.
    updates = [(classifier.W, classifier.W - learning_rate * g_W),
               (classifier.b, classifier.b - learning_rate * g_b)]

    # 这里定义一个训练模型的函数，这个函数是用来反复调用训练模型的，每一次调用函数，都会根据输入
    # 数据集计算这个模型当前参数下的损失函数，然后根据这个损失函数，用updates方式，对模型参数
    # 进行更新，第二次调用这个function，又会执行cost计算和updates计算，也就是通过反复调用
    # train_model，updates使得模型参数越来越优，cost损失函数的值越来越小。
    train_model = theano.function(
        inputs=[index],
        outputs=cost,
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        }
    )

这部分是训练模型的函数，在上面所有代码的基础上，这段代码是很好理解的，基本都是上面代码的组装。
基础模型已经搭建好了，等着另外写一篇博客仔细讲训练过程。

会另外具体讲解

这部分是预测函数，即我们在训练好模型以后，用这个模型做预测，使用的是train数据集

def predict():
    # load the saved model
    classifier = pickle.load(open('best_model.pkl'))

    # compile a predictor function
    predict_model = theano.function(
        inputs=[classifier.input],
        outputs=classifier.y_pred)

    # We can test it on some examples from test test
    dataset='mnist.pkl.gz'
    datasets = load_data(dataset)
    test_set_x, test_set_y = datasets[2]
    test_set_x = test_set_x.get_value()

    predicted_values = predict_model(test_set_x[:10])
    print("Predicted values for the first 10 examples in test set:")
    print(predicted_values)

if __name__ == '__main__':
    sgd_optimization_mnist()