深度学习算法实践5---线性回归算法实现

在学习了基本的Theano的概念之后，我们可以将这些知识用来进行一些简单的应用。在这篇文章中，我们将实现一个简旱的逻辑回归算法。逻辑回归的数学模型还是有些小复杂的，但是其基本概念却是非常简单的，因此，我们在这里只讨论基本概念，而不太涉及复杂的数学理论。
我们可以通过一个简单的例子来说明什么是逻辑回归算法，假设在三维空间中，有一组待分类的点，同时有一系平面，代表这些点应该属于的类别，我们将通过这些点到代表类别的平面的距离，来判断点所属于的类别。也就是说，对于一个点来，我们找到其距离最近的平面，那么我们就说这个点属于这个类别。上面的讨论是在三维空间下进行的，如果推广到多维空间，那么我们这里的平面，就变成了超平面，但是概念是类似的。
将上述描述换成数学语言，我们假设输入向量为x，其维数为D，输出类别为Y，共有N个类别，对于上面的分类问题，如果我们把问题简化，就变为对二维平面上的点，以及一系列代表类别的线，求距离该点最近的直线的问题，而直线在二维情况下可以表求为：y=wx + b，其中w为权重，b为偏移量，如果将上式推广到高维空间，则权值将变为一个矩阵，偏移量将变为一个矩阵，可以表求为Wx + b，我们将权值矩阵和偏移量向量称为模型的参数集。
有了上面的定义之后，我们可以得到如下数学表示：



其中i代表第i个类别，N代表总类别数，P(Y=i | x, W, b)代表是i类别的概率，对于哪个类别对应的概率最大，那么这个输入向量就应该属于哪个类别。在数学上可以表示为：

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
)
self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)
self.y_pred = T.argmax(self.p_y_given_x, axis=1)

在上面的代码中，我们用到的softmax函数，T.dot函数用于求向量的点积，即每个维度相乘再相加，T.argmax为求最大值。
逻辑回归算法中，可以定义对数自然函数，然后利用最大似然算法来找出最优解，但是计算机实现中，通常定义损失函数，将优化的目标定为求损失函数的最小值，对地逻辑回归问题，我们可以定义对数似然函数为：



我们直接将损失函数定义为似然函数的负数，如下所示：



我们现在的任务就变成了求损失函数的最小值问题，式中的

代表模型中的参数。虽然我们在定义损失函数时，用的是概率对数的和，但是我们由于要使用迷你批处理方式，因此我可定义损失函数为这些样本的平均值，即先求出每个数据样本的损失函数，然后求出这些损失函数的平均值，作为这个批次的损失函数。代码如下所示：

return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])

式中y.shape[0]是这个样本批次的个数，在这里我们高为n。T.arange(y.shape[0])为一个向量，值为0, 1, ..., n-1。
T.log(self.p_y_given_x是一个矩阵，行对应于第几个样本，列对应于该样本属于该类别概率的对数值，我们假定有|D|个类别，如下所示：
1 2 ... |D|
0 * * ... *
1 * * ... *
2 * * ... *
... ..............................
n-1 * * ... *
上面的式子虽然比较复杂，但是其本质是前面定义了一个矩阵，后边的表达式规定了矩阵的行数和列数，用于指导求其平均数。
有了上面的做准备工作，我们就可以定义一个逻辑回归的工具类了，其实这个类还是非常简单的。

class LogisticRegression(object):
def __init__(self, input, n_in, n_out):
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
)
self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)
self.y_pred = T.argmax(self.p_y_given_x, axis=1)
self.params = [self.W, self.b]
self.input = input

def negative_log_likelihood(self, y):
return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])

def errors(self, 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)
)
if y.dtype.startswith('int'):
return T.mean(T.neq(self.y_pred, y))
else:
raise NotImplementedError()

在有了之前的理论准备之后，上面的代码其实不难理解。
下面来讲解一下如何使用这个工具类来实现逻辑回归分类。在这里我们会使用MNIST手写数字识别库，该库包含60000个手写的数字图片，分辨率为28*28，采用黑底白字方式，黑色代表0，白色代表255。这个问题就变为输入向量维数为28*28=784，分为0~9类，共有10类的分类问题。顺便说一句，大家看到，很多机器学习算法，都采用MNIST的手写数字识别库，这是因为大家的算法需要一个标准化的东西进行性能和准确性方面的比较，因此人们需要一个标准的问题，大家的算法才有可比性，所以有人形象的说，MNIST手写数字识别库，就像生物学中的大肠杆菌一样，也是一个标准模型。
我们在使用逻辑回归算法之前，需要先确定一些参数，首先最重要的是学习率，学习率太小，模型收敛慢，而学习率太大，则不容易达到最优解，确定学习率是很多算法质量高低的重要方面，当将所有样本训练一遍，我们称其为一个epoch，我们不希望模型永远运行下去，如果始终达不到最优解，我们希望在到一定次数之后就停止执行，我们用epches来规定这个数值，由于我们是迷你批次方法来训练网络，所以我们需要批次所包含的样本量，在下面的代码中，所有这些参数均给出了缺省值。

def sgd_optimization_mnist(learning_rate=0.13, n_epochs=1000,
dataset='mnist.pkl.gz',
batch_size=600):
# 装入MNIST数据
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]
# 定义训练批次、验证批次、测试批次
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

index = T.lscalar()  # index to a [mini]batch
# 输入、输出向量定义
x = T.matrix('x')  # data, presented as rasterized images
y = T.ivector('y')  # labels, presented as 1D vector of [int] labels
# 分类器定义
classifier = LogisticRegression(input=x, n_in=28 * 28, n_out=10)
# 损失函数，学习算法使该值最小
cost = classifier.negative_log_likelihood(y)

# 测试模型定义
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]
}
)
# 梯度计算
g_W = T.grad(cost=cost, wrt=classifier.W)
g_b = T.grad(cost=cost, wrt=classifier.b)
# 更新参数值
updates = [(classifier.W, classifier.W - learning_rate * g_W),
(classifier.b, classifier.b - learning_rate * g_b)]
# 训练模型定义
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]
}
)
# 训练模型
patience = 5000
patience_increase = 2
improvement_threshold = 0.995
validation_frequency = min(n_train_batches, patience // 2)
best_validation_loss = numpy.inf
test_score = 0.
start_time = timeit.default_timer()
done_looping = False
epoch = 0
while (epoch < n_epochs) and (not done_looping):
epoch = epoch + 1
for minibatch_index in range(n_train_batches):
minibatch_avg_cost = train_model(minibatch_index)
iter = (epoch - 1) * n_train_batches + minibatch_index
if (iter + 1) % validation_frequency == 0:
validation_losses = [validate_model(i)
for i in range(n_valid_batches)]
this_validation_loss = numpy.mean(validation_losses)
if this_validation_loss < best_validation_loss:
if this_validation_loss < best_validation_loss *  \
improvement_threshold:
patience = max(patience, iter * patience_increase)
best_validation_loss = this_validation_loss
test_losses = [test_model(i)
for i in range(n_test_batches)]
test_score = numpy.mean(test_losses)
with open('best_model.pkl', 'wb') as f:
pickle.dump(classifier, f)
if patience <= iter:
done_looping = True
break
end_time = timeit.default_timer()

上面代码中，首先装入MNIST样本集，并划分为训练、验证、测试样本集，并形成培训的样本批次。接下来指定分类器为我们前面所定义的逻辑回归类，代价函数为负的对数似然函数。接下来定义了验证、测试、训练模型。接下来就是模型训练部分了。
在每个训练周期中（即epoch），依次循环每个批次，对于任一个批次，先对模型进行训练，然后看是否需要验证一下模型，如果需要则从验证模型中，求出代价函数的平均值，看看与记录的最小值是不是有显著改进，如果有则替换原来的最小值，并更新提早结束条件，并且将该阶段参数保存起来。
上面代码虽然在细节处理方面较为复杂，但是对于多数机器学习项目而言，基本是相同的，因此对于初学者而言，可以将其视为黑盒子。
其实对于实际应用而言，最重要的应该是对应用进行建模，然后装入数据对网络进行训练，所以下面重点讲述一下这部分相关的代码。在上面的代码中，调用了load_data函数，加载数据样本功能就是在这个函数中实现的。

def shared_dataset(data_xy, borrow=True):
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)
return shared_x, T.cast(shared_y, 'int32')

def load_data(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)
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)
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

有了上面这些代码之后，我们可以定义主函数，来训练这个逻辑回归算法了。

if __name__ == '__main__':
sgd_optimization_mnist()

运行这个程序，大约一两分钟之后，就可以得到一个训练好的逻辑回归算法模型了。
其实上面的代码完全是Github上面的一个开源的项目，大家完全可以参考这个项目。然而我们学习逻辑回归算法，并不仅仅是实现这一模型，利用MNIST数据集来进行训练，而是要解决我们的实际问题。但是网上绝大多数教程却仅止步于此，使得深度学习算法的学习入门非常困难。
所以在下一篇博文中，我们将以这个模型为基础，更换数据模型，来解决我们的实际问题。