机器学习-逻辑回归识别手写数字

写在前面

本博客是作者在学习机器学习基础时写下的总结，学习的资源为网易云课堂上吴恩达的机器学习课程，参考资料主要为stanford cs229课程的英文讲义,有兴趣的读者可以去网上下载原版的英文讲义来看。此外，由于本人为初学者，对知识点理解有限，因此文中有任何错误非常欢迎大家指出。最后，本文涉及的所有代码的完整版均会上传到github,欢迎大家交流。

本文主要内容课程的第三次作业的总结，由于课程作业的形式为填空的形式，很多实现细节容易遗漏，因此作者将作业中的octave/matlab程序用python进行重写，从而能够有更加深刻的理解。此外，本文在优化成本函数时采用了最简单的梯度下降法和scipy库中的minimize函数，并对minimiz函数的使用进行了简要说明。下面是正文：

1. 数据集

数据集是以matlab格式保存的5000张手写数字灰度图像，每张图像的大小为20*20像素，利用python可以打开.mat文件并转换为Numpy的数组，具体方法可参考https://blog.csdn.net/Bill_zhang5/article/details/79095985。下面是数据导入的部分代码与可视化后的数据集：

import scipy.io as scio
...
def dataRead():
    data = scio.loadmat('data.mat')
    X = data['X']
    y = data['y']
    return X, y

这里将5000张20*20像素图像组合成了一张1000*2000像素的图像进行显示。可以看到，0-9十个数字各有500张。

2. 逻辑回归

本文利用逻辑回归进行分类，输出的结果有10个类别。实际上，分10个类别和分2个类别的思想是一样的，也可以将分10个类别看做10次二分类。具体来说就是，利用10个假设函数分别对数据进行分类，则第i个假设函数则用于区分是否是第i个数字。下面进行具体分析：

假设函数

由于需要分出10个类别，因此需要有10个形式相同的假设函数，它们的不同之处在于回归之后的参数$\theta$不同。当某张图为数字i时，第i个假设函数的值应取到最大。

假设函数的形式为：
$$h_{\theta }(x)=g({\theta }^{T}x)$$
其中, $g(z)=\frac{1}{1+e^{-{\theta }^{T}x}}$为sigmoid函数。
输入为：$x=(1,x_1,x_2,...x_{400})$,
参数为：$\theta =({\theta }_0,{\theta }_1,...{\theta }_{400})$

注意：由于有10个假设函数，因此参数
$$\Theta =\left[\begin{array}{c}{\theta }^{(0)}\\ {\theta }^{(1)}\\ ...\\ {\theta }^{(9)}\end{array}\right]={\left[\begin{array}{c}{\theta }_0^{(0)},{\theta }_1^{(0)},...{\theta }_{400}^{(0)}\\ {\theta }_0^{(1)},{\theta }_1^{(1)},...{\theta }_{400}^{(1)}\\ ...\\ {\theta }_0^{(9)},{\theta }_1^{(9)},...{\theta }_{400}^{(9)}\end{array}\right]}_{10\times 401}$$
也有假设函数
$$H=\left[\begin{array}{c}{h}_{{\theta }^{(0)}}(x)\\ h_{{\theta }^{(1)}}(x)\\ ...\\ h_{{\theta }^{(9)}}(x)\end{array}\right]$$

成本函数

成本函数如下：
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\log (h_{\theta }(x^{(i)}))-(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))\right]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$
其中，$\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$为正则化项，注意j的取值为1到n而不是0到n。

梯度

对参数$\theta$求偏导得：
$$\frac{\partial J(\theta )}{\partial {\theta }_0}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_0^{(i)}$$
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\left[\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}\right]+\frac{\lambda }{m}{\theta }_{j}forj>0$$

矢量化

从上面的几个公式可以看出，在计算成本函数和梯度的时候有许多的乘法和求和运算。实际上，我们可以将上述的乘法和求和运算全部转换为矩阵运算，然后利用Numpy中高度优化的矩阵运算函数可以大大提高计算的效率。

在笔者看来，矢量化一个很重要的思想是把乘法+求和运算转换为类似于向量内积的形式，考虑下面公式：
$${\alpha }^T\beta ={\beta }^T\alpha =a^{(1)}{b}^{(1)}+a^{(2)}{b}^{(2)}+...+a^{(n)}{b}^{(n)}=\sum _{i=1}^n{a}^{(i)}{b}^{(i)}$$

其中，$\alpha =(a^{(1)},a^{(2)},...a^{(n)}{)}^T$,$\beta =(b^{(1)},b^{(2)},...b^{(n)}{)}^T$均为n维列向量。

类似的，对于成本函数：
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\log (h_{\theta }(x^{(i)}))-(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))\right]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$
可以通过矢量化转换为：
$$J(\theta )=-\frac{1}{m}\left[Y^T\log (h(X))-(1-Y{)}^T\log (1-h(X))\right]+\frac{\lambda }{2m}{\Theta }^T\Theta$$
注意：上式中正则化项不应取所有$\theta$，为了表述方便式子中就不加以区分了。对应matlab代码如下：

J = (-1/m) * (y' * log(sigmoid(X*theta)) + ...
              (1-y)' * log(1-sigmoid(X*theta))) ...
           + (lambda/(2*m)) * (theta(2:end)' * theta(2:end))

等价的python代码如下：

h = sigmoid(data.dot(theta))
J = (-1/self.m) * (label.T.dot(np.log(h)) +\
                   (1-label).T.dot(np.log(1-h)))\
                + (lam/2/m) * theta[1:].T.dot(theta[1:])   # regularization

同样的，对于矢量化之后的梯度为：
$$\frac{\partial J(\theta )}{\partial \theta }=\left[\frac{1}{m}(h(X)-Y{)}^{T}X\right]+\frac{\lambda }{m}\Theta$$
对应的matlab代码为：

grad = (1/m) * X' * (sigmoid(X * theta) - y);
grad(2:end) = grad(2:end) + (lambda/m) * theta(2:end); %regular term of gradience

python代码为：

h = self.sigmoid(data.dot(theta))
grad = (1 / self.m) * data.T.dot(h - label)
grad[1:] = grad[1:] + (self.lam / self.m * theta[1:])  # regularization

注：${\alpha }^T\beta ={\beta }^T\alpha$

利用梯度下降法的完整代码如下：

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


class MultiClassification(object):
    def __init__(self, data, label, classes, alpha, lam, iter):
        self.m = label.shape[0]	# size of trainint set	, 5000
        self.n = data.shape[1] + 1	# dimension of input, 401
        self.data = np.c_[np.ones((self.m, 1)), data]	# input data (image)
        self.label = label  # output , 0->9
        self.alpha = alpha	  # learning rate
        self.lam = lam		# regularization coefficient
        self.iter = iter	# iterations
        self.classes = classes	# number of classes (10)
        self.theta = np.zeros((classes, self.n))	# initial theta , 10 * 401

    def sigmoid(self, x):
        return 1 / (1 + np.exp(-x))  #sigmoid function

    def costFunc(self, i):
    	# argument i mean the ith 
        theta = self.theta[i, :].reshape(self.theta.shape[1], 1)
        label = (self.label == i + 1).astype('int16')  #(1)
        data = self.data
		# hypothesis function
        h = self.sigmoid(data.dot(theta))
        # regularized cost function
        J = (-1/self.m) * (label.T.dot(np.log(h))\
                            +(1-label).T.dot(np.log(1-h)))\
                     + (self.lam/2/self.m) * theta.T.dot(theta)   # regularization
        # hypothesis function
        grad = (1 / self.m) * data.T.dot(h - label)
        grad[1:] = grad[1:] + (self.lam / self.m * theta[1:])  # regularization
        return J, grad

    def regression(self):
        for i in range(self.classes):
            for j in range(self.iter):
                [J, grad] = self.costFunc(i)
                self.theta[i, :] = self.theta[i, :] - self.alpha * grad.T  # gradient descent
            print(J)
        return self.theta

    def predict(self):
    	# calc 10 Hypothesis functions and select the max one
    	# use argmax(1) to get the index of max_val of each row
        pred = self.sigmoid(self.data.dot(self.theta.T)).argmax(1) + 1
        # calc the accuracy of predict
        # the return value type of operator '==' is bool
        # use astype() to change bool into int16
        accuracy = np.mean(
            (self.label.flatten() == pred).astype('int16')) * 100
        print('accuracy:', accuracy)
        pass


def dataRead():
    data = scio.loadmat('data.mat')
    X = data['X']
    y = data['y']
    print(X.shape, y.shape)
    return X, y


def dataVisualize(data):  # combine them into a 1000*2000 image
    temp = data[0, :].reshape(20, 20)
    for i in range(4999):
        temp = np.c_[temp, data[i + 1, :].reshape(20, 20)]
    temp1 = temp[:, :2000]
    for i in range(49):
        temp1 = np.r_[temp1, temp[:, 2000 * (i + 1):2000 * (i + 2)]]
    plt.imshow(temp1)
    plt.show()


if __name__ == "__main__":
    [data, label] = dataRead()
    # dataVsualize(data)
    f = MultiClassification(data, label, 10, 1, 0.1, 2000)  # new a object
    f.regression()  # logistic regression
    f.predict()  # prediction

由于代码中有较详细的注释，因此这里不用进行过多解释，重要的有一下几点：

1. 输出Y的处理：label = (self.label == i + 1).astype('int16'),这里如果label=[1 2 3 4 3 2 5 5 1 5 1 5 3 1 1 3 9 8 7 9 6 6 6 9 1 3 6 5]且i=4，则该语句的结果为label = [0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]，然后进行二元分类即可。
2. 准确率的计算：accuracy = np.mean((self.label.flatten() == pred).astype('int16')) * 100直接相等的数据的均值，然后乘以100即可得到准确率。

下面是程序输出的结果：

[[ 0.01644014]]
[[ 0.05466642]]
[[ 0.0610377]]
[[ 0.03709997]]
[[ 0.05838671]]
[[ 0.02233645]]
[[ 0.03429286]]
[[ 0.08161475]]
[[ 0.07443368]]
[[ 0.01069821]]
accuracy: 95.32

通过10000次的迭代之后，计算得出准确率为95.32%。

虽然上述程序将所有计算公式进行了矢量化，但其使用的梯度下降法实际上是一种效率并不高的优化算法，因此上述程序迭代10000次在笔者的laptop上需要跑超过2分钟。下面介绍利用scipy中的minimize()函数对我们的成本函数进行优化。

3. 高级优化算法

scipy库这里就不多说了，想必大家都听说过，下面将用到scipy库中optimize模块中的minimize()函数对我们逻辑回归中的成本函数进行优化。这里是minimize()的官方介绍：https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.minimize.html

minmize()函数实现了课程中NG老师讲到的CG算法、BFGS算法、L-BFGS算法和很多其他高效的算法。函数原型为：
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)[source]

这里将笔者遇到的几个坑列出来进行总结和供大家参考：
1.需要分别写cost function和gradient function，它们的参数必须相同，它们分别需要传递给形参fun和jac。

2.cost function 返回值为成本函数的值，这里为标量；gradient function返回值为各参数的梯度，这里形状为(n,)而不是(n,1)

3.x0为需要优化的参数的初始值，这里是参数$\theta$的初始值，形状为(n, )

4.cost function和gradient functin的第一个参数应该用于接收上面传入的x0，也就是参数$\theta$。

5.args为需要额外传入的参数列表，可以用来传入训练数据什么的。注意args也会被传入cost function和gradient function

 def costFunc(theta, *args):
 	pass
 	return J
 def gradient(theta, *args):
	pass
	return grad.flatten()
  ...
 so.minimize(fun=costFunc
            ,x0=theta\
            , args=(args)\
            , method='TNC'\
            ,jac=gradient)['x'] 

6.minimize()函数的返回值不是优化后的参数$\theta$，而是一个详细信息的字典，需要我们从字典中提取出优化后的参数。

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as scio
import scipy.optimize as so


class MultiClassification(object):
    def __init__(self, data, label, classes, alpha, lam, iter):
        self.m = label.shape[0]
        self.n = data.shape[1] + 1
        self.data = np.c_[np.ones((self.m, 1)), data]
        self.label = label
        self.alpha = alpha
        self.lam = lam
        self.iter = iter
        self.classes = classes
        self.theta = np.zeros((classes, self.n))

    def sigmoid(self, x):
        return 1 / (1 + np.exp(-x))

    def costFunc(self, theta, i):
        theta = theta.reshape(theta.size, 1)
        label = (self.label == i+1).astype('int16')
        data = self.data
        h = self.sigmoid(data.dot(theta))
        J = (-1/self.m) * (label.T.dot(np.log(h)) +\
                          (1-label).T.dot(np.log(1-h)))\
                     + (self.lam/2/self.m) * theta[1:].T.dot(theta[1:])   # regularization
        print(J)
        return J

    def gradient(self, theta, i):
        theta = theta.reshape(theta.size, 1)
        label = (self.label == i+1).astype('int16')
        data = self.data
        
        h = self.sigmoid(data.dot(theta))
        grad = (1 / self.m) * data.T.dot(h - label)
        grad[1:] = grad[1:] + (self.lam / self.m * theta[1:])  # regularization
        # print(grad)
        return grad.flatten()

    def predict(self):
        pred = self.sigmoid(self.data.dot(self.theta.T)).argmax(1) + 1
        accuracy = np.mean(
            (self.label.flatten() == pred).astype('int16')) * 100
        print('accuracy:', accuracy)

    def regression(self):
        for i in range(self.classes):
            self.theta[i, :] = so.minimize(fun=self.costFunc, x0=self.theta[i, :]\
                                            , args=(i)\
                                            , method='TNC'\
                                            , jac = self.gradient)['x']

        return self.theta

   
def dataRead():
    data = scio.loadmat('data.mat')
    X = data['X']
    y = data['y']
    return X, y


def dataVisualize(data):  # combine them into a 1000*2000 image
    temp = data[0, :].reshape(20, 20)
    for i in range(4999):
        temp = np.c_[temp, data[i + 1, :].reshape(20, 20)]
    temp1 = temp[:, :2000]
    for i in range(49):
        temp1 = np.r_[temp1, temp[:, 2000 * (i + 1):2000 * (i + 2)]]
    plt.imshow(temp1, cmap='gray')
    plt.show()


if __name__ == "__main__":
    [data, label] = dataRead()
    dataVisualize(data)
    f = MultiClassification(data, label, 10, 1, 0.1, 2000)
    f.regression()
    f.predict()

运行结果为：

accuracy: 96.46

可以看到，模型在训练集上的准确率达96.46%，且该程序在笔者的laptop上运行时间不超过20秒，比之前用梯度下降法快了很多倍。因此可以看出利用高性能的优化算法可以让训练速度大大提高且不需要手动调试学习率$\alpha$等参数。

值得注意的是，上面optimize()函数中选择的是’TNC’算法，它的全称是"truncated Newton algorithm"，中文翻译是截断牛顿算法。经过笔者测试，该算法在本文的模型下比CG算法、BFGS算法、L-BFGS算法效果更好，效率更高。这三种算法中，CG(共轭梯度法)算法比其他两种算法性能要好。(注：仅在本模型下）

另外，笔者作为一个初学者只能看到上述结果的表象，如果存在错误欢迎指出！

最后，本文的所有代码和数据集都在前言中的github中。