梯度下降Python实现

梯度下降

• 公式
  
  梯度下降就是根据当下所在的位置寻找一个最陡峭的地方往下，在这一点的切线上我认为有两个方向，一个向下一个向上，把向下的方向分为水平方向的向量和垂直向下的向量，水平方向的向量的长度就是要跨越的距离，与η有关（η不宜太大，不然会错过最陡峭的地方，也不宜太小），水平跨越后的结果就是上面的公式。（梯度上升就是把向上的方向分解）
  

（图画的稍微潦草看懂即可）如何求最陡峭的地方就要用到数学上的微分了。
1、单微分：

2、多微分：

梯度一般来说就是有多变量的微分函数，在这一点上也可以分为两个方向的向量，向量有方向，梯度的方向就指出了函数在给定点的上升最快的方向。

代码及公式推导（参考自网上博客）

下面我们将用python实现一个简单的梯度下降算法。

我们将用梯度下降法来拟合出这条直线！

首先，我们需要定义一个代价函数，在此我们选用均方误差代价函数

此公示中

• m是数据集中点的个数
• 1/2是一个常量，这样是为了在求梯度的时候，二次方乘下来就和这里的1/2抵消了，自然就没有多余的常数系数，方便后续的计算，同时对结果不会有影响
• y 是数据集中每个点的真实y坐标的值
• h 是我们的预测函数，根据每一个输入x，根据Θ 计算得到预测的y值，即

我们可以根据代价函数看到，代价函数中的变量有两个，所以是一个多变量的梯度下降问题，求解出代价函数的梯度，也就是分别对两个变量进行微分

明确了代价函数和梯度，以及预测的函数形式。我们就可以开始编写代码了。但在这之前，需要说明一点，就是为了方便代码的编写，我们会将所有的公式都转换为矩阵的形式，python中计算矩阵是非常方便的，同时代码也会变得非常的简洁。

为了转换为矩阵的计算，我们观察到预测函数的形式

我们有两个变量，为了对这个公式进行矩阵化，我们可以给每一个点x增加一维，这一维的值固定为1，这一维将会乘到Θ0上。这样就方便我们统一矩阵化的计算

然后我们将代价函数和梯度转化为矩阵向量相乘的形式
接下来编写代码

#首先，我们需要定义数据集和学习率
from numpy import *
import matplotlib.pyplot as plt
# 数据集大小 即20个数据点
m = 20
# x的坐标以及对应的矩阵
X0 = ones((m, 1))  # 生成一个m行1列的向量，也就是x0，全是1
X1 = arange(1, m+1).reshape(m, 1)  # 生成一个m行1列的向量，也就是x1，从1到m
X = hstack((X0, X1))  # 按照列堆叠形成数组，其实就是样本数据
# 对应的y坐标
Y = array([
    3, 4, 5, 5, 2, 4, 7, 8, 11, 8, 12,
    11, 13, 13, 16, 17, 18, 17, 19, 21
]).reshape(m, 1)
# 步长
alpha = 0.01

#接下来我们以矩阵向量的形式定义代价函数和代价函数的梯度
# 定义代价函数
def cost_function(theta, X, Y):
    diff = dot(X, theta) - Y  # dot() 数组需要像矩阵那样相乘，就需要用到dot()
    return (1/(2*m)) * dot(diff.transpose(), diff)
# 定义代价函数对应的梯度函数
def gradient_function(theta, X, Y):
    diff = dot(X, theta) - Y
    return (1/m) * dot(X.transpose(), diff)

#最后就是算法的核心部分，梯度下降迭代计算
# 梯度下降迭代
def gradient_descent(X, Y, alpha):
    theta = array([1, 1]).reshape(2, 1)
    gradient = gradient_function(theta, X, Y)
    while not all(abs(gradient) <= 1e-5):
        theta = theta - alpha * gradient
        gradient = gradient_function(theta, X, Y)
    return theta
optimal = gradient_descent(X, Y, alpha)
print('optimal:', optimal)
print('cost function:', cost_function(optimal, X, Y)[0][0])

当梯度小于1e-5时，说明已经进入了比较平滑的状态，类似于山谷的状态，这时候再继续迭代效果也不大了，所以这个时候可以退出循环！
运行代码，计算得到的结果如下：

print('optimal:', optimal)
print('cost function:', cost_function(optimal, X, Y)[0][0])

# 根据数据画出对应的图像
def plot(X, Y, theta):
    ax = plt.subplot(111)
    ax.scatter(X, Y, s=20, c="green", marker="s")
    plt.xlabel("X")
    plt.ylabel("Y")
    x = arange(0, 22, 0.01)  # x的范围
    y = theta[0] + theta[1]*x
    ax.plot(x, y)
    plt.show()
plot(X1, Y, optimal)

运行结果