机器学习实战之Logistic回归——最优化算法的讨论

写这篇文章的时候，我还蛮想借鉴《机器学习实战》一书中这一章节的开头语：这会是激动人心的一章，因为我们将首次接触到最优化算法。说得真好啊，所以希望各位读者在阅读这篇文章时能按耐住自己的激动之情好好看下去(哈哈哈哈，有点欠揍)。言归正传，这一章节的核心就是讨论Logistic回归中使用到的一些最优化算法，比如梯度下降，牛顿法和拟牛顿法，但是书中所给的基本上都是用梯度上升算法实现的案例。正所谓不走寻常路，我将在下面好好讨论梯度上升，梯度下降算法，牛顿法和拟牛顿法，最后还将会借鉴神经网络中的优化算法来给各位读者做延伸阅读。因为我还是那句老话：不敷衍，不忽悠。只是希望各位能有些感触(嘿嘿，是不是觉得我形象立马高大了起来)

1. 梯度下降相关概念简介

其实我是最烦写概念的了(其实还是肚子里没墨水)，所以这一段概念简介我基本上就是借鉴了我认为写的很好的一位博主的文章，大家很容易看懂。https://www.jianshu.com/p/c7e642877b0e这个是这篇文章的链接，大家如果想好好弄懂相关概念可以去这篇文章看看(但不要看完那篇文章之后就抛弃我了啊，好好看后面的总结不会让你失望的)

梯度下降的场景假设

梯度下降法的基本思想可以类比为一个下山的过程。假设这样一个场景：一个人被困在山上，需要从山上下来(i.e. 找到山的最低点，也就是山谷)。但此时山上的浓雾很大，导致可视度很低。因此，下山的路径就无法确定，他必须利用自己周围的信息去找到下山的路径。这个时候，他就可以利用梯度下降算法来帮助自己下山。具体来说就是，以他当前的所处的位置为基准，寻找这个位置最陡峭的地方，然后朝着山的高度下降的地方走，同理，如果我们的目标是上山，也就是爬到山顶，那么此时应该是朝着最陡峭的方向往上走。然后每走一段距离，都反复采用同一个方法，最后就能成功的抵达山谷。

我们同时可以假设这座山最陡峭的地方是无法通过肉眼立马观察出来的，而是需要一个复杂的工具来测量，同时，这个人此时正好拥有测量出最陡峭方向的能力。所以，此人每走一段距离，都需要一段时间来测量所在位置最陡峭的方向，这是比较耗时的。那么为了在太阳下山之前到达山底，就要尽可能的减少测量方向的次数。这是一个两难的选择，如果测量的频繁，可以保证下山的方向是绝对正确的，但又非常耗时，如果测量的过少，又有偏离轨道的风险。所以需要找到一个合适的测量方向的频率，来确保下山的方向不错误，同时又不至于耗时太多！

梯度下降

梯度下降的基本过程就和下山的场景很类似。

首先，我们有一个可微分的函数。这个函数就代表着一座山。我们的目标就是找到这个函数的最小值，也就是山底。根据之前的场景假设，最快的下山的方式就是找到当前位置最陡峭的方向，然后沿着此方向向下走，对应到函数中，就是找到给定点的梯度，然后朝着梯度相反的方向，就能让函数值下降的最快！因为梯度的方向就是函数之变化最快的方向(在后面会详细解释)
所以，我们重复利用这个方法，反复求取梯度，最后就能到达局部的最小值，这就类似于我们下山的过程。而求取梯度就确定了最陡峭的方向，也就是场景中测量方向的手段

梯度

梯度实际上就是多变量微分的一般化。
下面这个例子：

我们可以看到，梯度就是分别对每个变量进行微分，然后用逗号分割开，梯度是用<>包括起来，说明梯度其实一个向量。

梯度是微积分中一个很重要的概念，之前提到过梯度的意义

在单变量的函数中，梯度其实就是函数的微分，代表着函数在某个给定点的切线的斜率
在多变量函数中，梯度是一个向量，向量有方向，梯度的方向就指出了函数在给定点的上升最快的方向

这也就说明了为什么我们需要千方百计的求取梯度！我们需要到达山底，就需要在每一步观测到此时最陡峭的地方，梯度就恰巧告诉了我们这个方向。梯度的方向是函数在给定点上升最快的方向，那么梯度的反方向就是函数在给定点下降最快的方向，这正是我们所需要的。所以我们只要沿着梯度的方向一直走，就能走到局部的最低点！

梯度下降算法的数学解释

上面我们花了大量的篇幅介绍梯度下降算法的基本思想和场景假设，以及梯度的概念和思想。下面我们就开始从数学上解释梯度下降算法的计算过程和思想！

此公式的意义是：J是关于Θ的一个函数，我们当前所处的位置为Θ0点，要从这个点走到J的最小值点，也就是山底。首先我们先确定前进的方向，也就是梯度的反向，然后走一段距离的步长，也就是α，走完这个段步长，就到达了Θ1这个点！

下面就这个公式的几个常见的疑问：

α是什么含义？
α在梯度下降算法中被称作为学习率或者步长，意味着我们可以通过α来控制每一步走的距离，以保证不要步子跨的太大扯着蛋，哈哈，其实就是不要走太快，错过了最低点。同时也要保证不要走的太慢，导致太阳下山了，还没有走到山下。所以α的选择在梯度下降法中往往是很重要的！α不能太大也不能太小，太小的话，可能导致迟迟走不到最低点，太大的话，会导致错过最低点！

为什么要梯度要乘以一个负号？
梯度前加一个负号，就意味着朝着梯度相反的方向前进！我们在前文提到，梯度的方向实际就是函数在此点上升最快的方向！而我们需要朝着下降最快的方向走，自然就是负的梯度的方向，所以此处需要加上负号

2. 谈论梯度上升和梯度下降

这么吧，我觉得我还是先把梯度上升和梯度下降的两个函数分别贴出来，先让大家有个比较形象的比较，我再来讲讲他们之间的区别(其实有个屁的区别)

# 梯度上升优化函数
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = np.ones((n, 1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

# 梯度下降优化函数
def gradient_descent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    weights = np.ones((n, 1))
    error = sigmoid(dataMatrix * weights) - labelMat
    gradient = dataMatrix.transpose() * error/m
    while np.all(np.absolute(gradient) >= 1e-8):
        h = sigmoid(dataMatrix * weights)
        weights -= alpha * gradient
        error = h - labelMat
        gradient = dataMatrix.transpose() * error/m
    return weights

下面这张图是我在我的书上做的笔记(本人较懒，相信你们早就看出来了，哎) 关于梯度▽f(x, y)的具体求导公式大家可以看看周志华老师的西瓜书和李航老师的《统计学习方法》(这两本书个人认为是我们国内机器学习的圣经了) 里面有详细的推导公式，这儿我就直接用了，希望大家理解

看完上面两个函数之后，相信大家会有跟我一样的感受，这两个算法就没啥区别，只不过是里面加减法变了一些。这时候有人就会站出来说：weights权重的变化一个是加法(梯度上升)即：weights += alphs*gradient ; 一个是减法(梯度下降)即：weights -= alphs*gradient 但是只有大家细心一点就会发现再计算error(即预测值与真实值之间的差值)时顺序是相反的，在梯度上升算法中error = labelMat - h 即是用真实值-预测值；而在梯度下降算法中是用error = h - labelMat 即是用预测值-真实值。了解到这点后大家是不是豁然开朗了，并且想爆句粗口：sha b... TM欺骗我感情最后再用书上的一句原话结束这一小节

梯度上升算法一般用来求最大化似然函数的最大值；梯度下降算法一般用来求最小化损失函数的最小值

3. 浅析随机梯度下降

个人觉得：扯了这么多关于梯度下降的内容，不进一步扯扯随机梯度下降的东西有点太可惜了。话不多说看下文

我们知道梯度上升法每次更新回归系数都需要遍历整个数据集，当样本数量较小时，该方法尚可，但是当样本数据集非常大且特征非常多时，那么随机梯度下降法的计算复杂度就会特别高(不信，你可以拿个十几兆的数据集在普通笔记本上跑，会让你崩溃的)。一种改进的方法是一次仅用一个样本点来更新回归系数，即岁集梯度上升法。由于可以在新样本到来时对分类器进行增量式更新，因此随机梯度上升法是一个在线学习算法。

随机梯度上升法可以写成如下伪代码：

所有回归系数初始化为1

对数据集每个样本

　　 计算该样本的梯度

　　 使用alpha*gradient更新回顾系数值

返回回归系数值

stocGradAscent0()函数如下：

# 随机梯度上升算法
def stocGradAscent0(dataMatrix, classLabels):
    m, n = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

我们知道，评判一个优化算法的优劣的可靠方法是看其是否收敛，也就是说参数的值是否达到稳定值。此外，当参数值接近稳定时，仍然可能会出现一些小的周期性的波动。这种情况发生的原因是样本集中存在一些不能正确分类的样本点(数据集并非线性可分)，所以这些点在每次迭代时会引发系数的剧烈改变，造成周期性的波动。显然我们希望算法能够避免来回波动，从而收敛到某个值，并且收敛速度也要足够快。

　　为此，需要对上述随机梯度上升法代码进行适当修改，代码如下：

# 改进后的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01
            randIndex = int(random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

上述代码中有两处改进的地方：

1 alpha在每次迭代更新是都会调整，这会缓解数据波动或者高频运动。此外，alpha还有一个常数项，目的是为了保证在多次迭代后仍然对新数据具有一定的影响，如果要处理的问题是动态变化的，可以适当加大该常数项，从而确保新的值获得更大的回归系数。

2 第二个改进的地方是选择随机的样本对参数进行更新，由于增加了随机性，这就防止参数值发生周期性的波动。

3 并且采用上述改进算法后，收敛速度更快

4.牛顿法

《机器学习实战》这本书上是没有介绍牛顿法以及后面的拟牛顿法，我是看到了李航老师写的《统计学习方法》上有介绍到了这两种方法，据说是效果不错，所以想来学习学习，也想分享给大家

4.1 牛顿法公式推导

看公式特别头疼，我先贴出源码，大家看看源码后，我再指出一些重点

# 计算海森矩阵
def computeHessianMatrix(dataMatIn):
    dataMatrix = np.mat(dataMatIn)
    hessianMatrix = dataMatrix.transpose().dot(dataMatrix)
    return hessianMatrix.I


# 牛顿法优化函数
def newtonMethod(dataMatIn, classLabels, numIter=100):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(dataMatrix)
    weights = np.ones((n, 1))
    alpha = 0.01
    hessianMatrix = computeHessianMatrix(dataMatrix)
    for k in range(numIter):
        h = sigmoid(dataMatrix * weights)
        error = labelMat - h
        weights -= (dataMatrix * hessianMatrix).transpose() * error
        # 加上了alpha学习率就变成了阻尼牛顿法
        # weights -= (dataMatrix * hessianMatrix).transpose() * error * alpha
    return weights

1. 不知道大家发现没有：牛顿法是要求二阶导数的，即海森矩阵。而我给出的海森矩阵：

hessianMatrix = dataMatrix.transpose().dot(dataMatrix)

特别简洁，至于为什么这么简洁。我目前还无法解释(没办法，水平太低，看了好几天都没看懂) 只能说大家记着就行

2. 前面说到了：梯度下降是求一阶导数的，牛顿法是求二阶导数的，所以牛顿法的时间复杂度是要低于梯度下降法的。现实也确实如此：一般的，梯度下降法迭代次数要到10000次左右才能有较好的分类效果，而牛顿法只需迭代到500次左右就能有和梯度下降法差不多的分类效果；还有至于牛顿法求二阶导数为什么会有这么好的效果，我理解是：一阶导数如果理解成梯度的话(即当前位置变化最快的方向) 那么二阶导数你可以理解为：把当前梯度当做参数，即求当前梯度变化最快的方向，相当于是从长远角度来考虑了。这样的话可以避免局部最优解

3. 我看了很多资料，虽然牛顿法有很好的效果，但是在神经网络中使用更多的却是梯度下降最优化方法。一是涉及参数较少，所以求解简单，因为牛顿法需要求海森矩阵，当参数较多时，海森矩阵也会变得异常庞大，效率反而降低。但是我觉得多了解这些方法对大家也是没坏处的，是吧

5. 神经网络延伸扩展

在第3节的最后部分我提到了改进后的随机梯度上升算法，你会发现：整到最后有太多太多的参数需要随着训练次数的增加而改变了，不然效果就没那么好了。比较懒的朋友(比如像我这类)就会觉得很烦，甚至想抛弃这种算法，由于接触了神经网络的相关内容，会发现使用Tensorflow工具来实现梯度下降(梯度上升)算法简直不要太easy了。都有现成的包，只需要一两行代码套公式就能轻轻松松实现了。

# 计算交叉熵的损失函数
    cross_entropy = tf.nn.sparse_softmax_cross_entropy_with_logits(logits=y, labels=tf.argmax(y_, 1))
    cross_entropy_mean = tf.reduce_mean(cross_entropy)
    loss = cross_entropy_mean + tf.add_n(tf.get_collection('losses'))
    learning_rate = tf.train.exponential_decay(LEARNING_RATE_BASE, global_step,
                                               mnist.train.num_examples / BATCH_SIZE, LEARNING_RATE_DECAY)
    train_step = tf.train.GradientDescentOptimizer(learning_rate).minimize(loss, global_step=global_step)

上面是我贴出来几行核心的代码，其实每行代码所要求的值很容易看懂。是不是觉得很easy呢，感兴趣的朋友也可以去阅读相关书籍，由于篇幅原因(其实还是那个原因，肚子里没墨水，哈哈哈哈哈哈)我在这就不累赘了，相当于抛砖引玉吧