1. 梯度下降法简介

1-1

以下是定义了一个损失函数以后，参数theta对应的损失函数J的值对应的示例图，我们需要找到使得损失函数值J取得最小值对应的theta（这里是二维平面，也就是我们的参数只有一个）

在直线方程中，导数代表斜率
在曲线方程中，导数代表切线斜率
导数代表theta单位变化时，J相应的变化

1-2
1-3

η太小，会减慢收敛学习速度

1-4
η太大，甚至导致不收敛
1-5

其他注意事项

• 并不是所有函数都有唯一的极值点

  
  1-6
  

  解决方案：

  • 多次运行，随机化初始点
  • 梯度下降法的初始点也是一个超参数

1-7

---

2 梯度下降法模拟

2.1 实现

import numpy as np
import matplotlib.pyplot as plt
# 简单模拟一个损失函数
plot_x = np.linspace(-1,6,141)
plot_y = (plot_x-2.5)**2-1
# 绘制我们模拟的损失函数
plt.plot(plot_x,plot_y)

2.1-1

def dJ(theta):
    """损失函数的导数"""
    return 2*(theta-2.5)

def J(theta):
    """损失函数"""
    try:
        return (theta-2.5)**2-1
    except:
        return float('inf')

def gradient_descent(initial_theta,eta,n_iters = 1e4,epsilon=1e-8):
    """
    梯度下降法封装
    initial_theta:初始化的theta值
    eta:学习率η
    n_iters: 最大循环次数
    epsilon: 精度
    """
    theta = initial_theta
    # theta_history 保存theta的变化值
    theta_history.append(initial_theta)
    i_iters = 0
    
    while i_iters<n_iters:
        """
        如果theta两次变化之间的损失函数值的变化小于我们定义的精度
        则可以说明我们已经找到了最低的损失函数值和对应的theta
        
        如果循环次数超过了我们设置的循环次数，
        则说明可能由于η设置的过大导致无止境的循环
        """
        gradient = dJ(theta)
        last_theta = theta
        theta = theta - eta * gradient
        theta_history.append(theta)
        
        if (abs(J(theta)-J(last_theta)) < epsilon):
            break
        i_iters += 1

def plot_theta_history():
    plt.plot(plot_x,J(plot_x))
    plt.plot(np.array(theta_history),J(np.array(theta_history)),color='r',marker='+')
    print("the size of theta_history is %d"%len(theta_history))

2.2 使用不同η学习率测试并观察我们的梯度下降法的结果

eta = 0.1
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-1

eta = 0.01
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-2

eta = 0.001
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-3

eta = 0.8
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-4

可以发现，只要η不超过一个限度，我们编写的函数都可以在有限次数之后找到最优解，并且η越小，学习的次数越多
下面来看一下，如果eta取较大值1.1，会出现什么情况

eta = 1.1
theta_history = []
gradient_descent(0.,eta)
# 数据量太大会报错
# plot_theta_history()
print(len(theta_history))
# 输出10001
theta_history[-1]
# 输出 nan（not a number）

可以看出当我们的eta取1.1，函数会循环直至终止，这是由于，我们的η设置过大，导致每次循环过后，损失函数j的值都向大的方向变化

eta = 1.1
theta_history = []
gradient_descent(0.,eta,n_iters = 10)
plot_theta_history()

2.2-5

---

3 多元线性回归中的梯度下降法

3-1

一个三维空间中的梯度下降法（x,y为系数，z为损失函数）

3-2

推导过程

3-3
3-4
上面推导出的式子的大小是和样本数有关的，m越大，结果越大，这是不合理的，我们希望和m无关
3-5

---

4 线性回归中的梯度下降法的实现

4.1 比较笨的方法实现

def fit_gd(self, X_train, y_train, eta=0.01, n_iters = 1e4):
        """根据训练数据集X_train,y_train, 使用梯度下降法训练Linear Regression 模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"

        def J(theta, X_b, y):
            try:
                return np.sum((y - X_b.dot(theta))**2) / len(X_b)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            res = np.empty(len(theta))
            res[0] = np.sum(X_b.dot(theta) - y)

            for i in range(1, len(theta)):
                res[i] = np.sum((X_b.dot(theta) - y).dot(X_b[:, i]))

            return res * 2 / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=n_iters, epsilon=1e-8):
            """
            梯度下降法封装
            X_b: X特征矩阵
            y: 结果向量
            initial_theta:初始化的theta值
            eta:学习率η
            n_iters: 最大循环次数
            epsilon: 精度
            """
            theta = initial_theta
            i_iters = 0

            while i_iters < n_iters:
                """
                如果theta两次变化之间的损失函数值的变化小于我们定义的精度
                则可以说明我们已经找到了最低的损失函数值和对应的theta
                
                如果循环次数超过了我们设置的循环次数，
                则说明可能由于η设置的过大导致无止境的循环
                """
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient

                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break

                i_iters += 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, initial_theta, eta)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

4.2 测试我们的算法

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x = 2*np.random.random(size = 100)
# 定义截距为4 斜率为3
y = x * 3. + 4. + np.random.normal(size=100)
X = x.reshape(-1,1)
plt.scatter(x,y)

4.2-1

from machine_learning.LinearRegression import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit_gd(X,y)

lin_reg.interception_
# 4.021457858204859
lin_reg.coef_
# array([3.00706277])

4.3 向量化

4.3-1
4.3-2
4.3-3

修改之前的求导函数

def dJ(theta, X_b, y):
            # res = np.empty(len(theta))
            # res[0] = np.sum(X_b.dot(theta) - y)
            #
            # for i in range(1, len(theta)):
            #     res[i] = np.sum((X_b.dot(theta) - y).dot(X_b[:, i]))
            #
            # return res * 2 / len(X_b)
            return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)

使用真实的数据测试

4.3-4

使用真实的数据，调整eta和iters，要么由于eta太小导致无法得出真实的结果，要么由于eta太大导致训练时间加长，这是由于数据的规模在不同的特征上不同，所以我们需要对数据进行归一化

4.4 数据归一化

4.4-1 4.4-2 4.4-3

如果样本数非常多，那么即使使用梯度下降法也会导致速度比较慢，因为在梯度下降法中，每一个样本都要参与运算。这时候需要采用随机梯度下降法，我们将在下一小节进行介绍

---

5 随机梯度下降法

5.1 随机梯度下降法介绍

5.1

批量梯度下降法带来的一个问题是η的值需要设置的比较小，在样本数比较多的时候导致不是速度特别慢，这时候观察随机梯度下降法损失函数的求导公式，可以发现，我们对每一个Xb都做了求和操作，又在最外面除以了m，那么可以考虑将求和和除以m的两个运算约掉，采用每次使用一个随机的Xb

5.2
5.3

由于我们使用的事随机梯度下降法，所以导致我们的最终结果不会像批量梯度下降法一样准确的朝着一个方向运算，而是曲线行下降，这时候我们就希望，越到下面，η值相应减小，事运算次数变多，从而精确计算结果

5-4
这里使用了模拟退火的思想
5-5

5.2 随机梯度下降法实现

def dJ_sgd(theta,X_b_i,y_i):
    return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2

def sgd(X_b,y,initial_theta,n_iters):
    
    t0 = 5
    t1 = 50
    
    def learning_rate(t):
        return t0 / (t + t1)
    
    theta = initial_theta
    for cur_iter in range(n_iters):
        rand_i = np.random.randint(len(X_b))
        gradient = dJ_sgd(theta,X_b[rand_i],y[rand_i])
        theta = theta - learning_rate(cur_iter) * gradient
    return theta

        %%time
        eta = 0.01
        X_b = np.hstack([np.ones((len(X), 1)), X])
        initial_theta = np.zeros(X_b.shape[1])
        # 随机的检查了3分之一个样本总量的样本
        _theta = sgd(X_b, y, initial_theta, n_iters=len(X_b)//3)
        # 输出 CPU times: user 318 ms, sys: 5.22 ms, total: 323 ms
        # Wall time: 337 ms
        # _theta: array([3.03182269, 3.93118623])

5.3 随机梯度下降法的封装和测试

    def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
        """
        根据训练数据集X_train, y_train, 使用随机梯度下降法训练Linear Regression模型
        :param X_train:
        :param y_train:
        :param n_iters: 在随机梯度下降法中，n_iters代表所有的样本会被看几圈
        :param t0:
        :param t1:
        :return:
        """
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        assert n_iters >= 1

        def dJ_sgd(theta, X_b_i, y_i):
            """
            去X_b,y 中的随机一个元素进行导数公式的计算
            :param theta:
            :param X_b_i:
            :param y_i:
            :return:
            """
            return X_b_i * (X_b_i.dot(theta) - y_i) * 2.

        def sgd(X_b, y, initial_theta, n_iters, t0=5, t1=50):

            def learning_rate(t):
                """
                计算学习率，t1 为了减慢变化速度，t0为了增加随机性
                :param t: 第t次循环
                :return:
                """
                return t0 / (t + t1)

            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                # 对X_b进行一个乱序的排序
                indexes = np.random.permutation(m)
                X_b_new = X_b[indexes]
                y_new = y[indexes]

                # 对整个数据集看一遍
                for i in range(m):
                    gradient = dJ_sgd(theta, X_b_new[i], y_new[i])
                    theta = theta - learning_rate(cur_iter * m + i) * gradient

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.random.randn(X_b.shape[1])
        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

模拟数据进行测试

5.3-1

真实数据波士顿房价进行测试

5.3-2

5.4 使用Sklearn中的 随机梯度下降法

5.4-1

需要注意的是sklearn中的梯度下降法比我们自己的算法要复杂的多，性能和计算准确度上都比我们的要好，我们的算法只是用来演示过程，具体生产上的使用还是应该使用Sklearn提供的

---

6 梯度下降法 的调试

6.1 梯度下降法调试的原理

可能我们计算出梯度下降法的公式，并使用python编程实现，预测的过程中程序并没有报错，但是可能我们需要求的梯度的结果是错误的，这个时候需要怎么样去调试发现错误呢。

首先以二维坐标平面为例，一个点（O）的导数就是曲线在这个点的切线的斜率，在这个点两侧各取一个点（AB），那么AB两点对应的直线的斜率应该大体等于O的切线的斜率，并且这A和B的距离越近，那么两条直线的斜率就越接近
事实上，这也正是导数的定义，当函数y=f(x)的自变量x在一点x₀上产生一个增量Δx时，函数输出值的增量Δy与自变量增量Δx的比值在Δx趋于0时的极限a如果存在，a即为在x₀处的导数，记作f'(x₀)或df(x₀)/dx

6-1

扩展到多维维度则如下

6-2

梯度下降法调试的实现

np.random.seed(666)
X = np.random.normal(size=(1000,10))
X_b = np.hstack([np.ones((len(X),1)),X])

# 真实的θ值
true_theta = np.arange(1,12,dtype=float)
# np.random.normal(size=1000) 添加噪音
y = X_b.dot(true_theta) + np.random.normal(size=1000)

# 实现J(θ)函数
def J(theta,X_b,y):
    try:
        return np.sum((y-X_b.dot(theta))**2)/len(X_b)
    except:
        return float('inf')
    
# 实现数学推导出的dJ(θ)
def d_J_main(theta,X_b,y):
    return X_b.T.dot(X_b.dot(theta) - y) * 2. /len(X_b)

# 实现debug模式的dJ(θ)
def d_J_debug(theta,X_b,y,epsilon=0.01):
    res = np.empty(len(theta))
    for i in range(len(theta)):
        theta_1 = theta.copy()
        theta_1[i] += epsilon
        theta_2 = theta.copy()
        theta_2[i] -= epsilon
        res[i] = (J(theta_1,X_b,y)-J(theta_2,X_b,y))/(2*epsilon)
    return res

# 批量梯度下降法，d_J为求导函数，作为一个参数传入，用于切换求导策略
def gradient_descent(d_J,X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
            theta = initial_theta
            i_iters = 0

            while i_iters < n_iters:
                gradient = d_J(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient

                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break

                i_iters += 1

            return theta

X_b = np.hstack([np.ones((len(X), 1)), X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01
# 使用d_J_debug调试模式求出theta
%time theta = gradient_descent(d_J_debug,X_b, y, initial_theta, eta)
print(theta)
# 使用数学解求出theta
%time theta = gradient_descent(d_J_main,X_b, y, initial_theta, eta)
print(theta)

# 输出结果
CPU times: user 531 ms, sys: 214 ms, total: 745 ms
Wall time: 613 ms
[ 0.94575233  1.98082712  3.06882065  3.94835863  4.97139932  5.9859077
  7.01077392  7.99250414  8.99151383  9.97525811 10.99758484]
CPU times: user 67 ms, sys: 27.6 ms, total: 94.6 ms
Wall time: 76.7 ms
[ 0.94575233  1.98082712  3.06882065  3.94835863  4.97139932  5.9859077
  7.01077392  7.99250414  8.99151383  9.97525811 10.99758484]

由此可以看出，我们的d_J_debug和d_J_main的结果是相近的，所以我们的d_J_main的数学推导是没问题的。
我们可以在真正的机器学习之前，先使用d_J_debug这种调试方式来验证一下我们的d_J_main的结果是否正确，然后再进行机器学习。

d_J_debug是通用的，可以放在任何求导的debug过程中，所以可以作为我们机器学习的工具箱来使用

---

7.梯度下降法的总结

7.1 小批量

• 批量梯度下降法
• 随机梯度下降法
  下面来看下二者的对比

维度	批量梯度下降法	随机梯度下降法
计算方式	每次对所有的样本看一遍才可以计算出梯度	每一次只需观察一个样本
速度	慢	快
稳定性	高，一定可以先向损失函数下降的方式前进	低，每一次的方式不确定，甚至向反方向前进

综合二者的优缺点，有一种新的梯度下降法

7-1

小批量梯度下降法：即，我们每一次不看全部样本那么多，也不是只看一次样本那么少，每次只看k个样本
对于小批量梯度下降法，由多了一个超参数

def fit_lit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50,k=10):
        """
        根据训练数据集X_train, y_train, 使用随机梯度下降法训练Linear Regression模型
        :param X_train:
        :param y_train:
        :param n_iters: 在随机梯度下降法中，n_iters代表所有的样本会被看几圈
        :param t0:
        :param t1:
        :param k: 小批量随机下降法的超参数k
        :return:
        """
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        assert n_iters >= 1

        def dJ_sgd(theta, X_b_k, y_k):
            """
            去X_b,y 中的随机选择k个元素进行导数公式的计算
            :param theta:
            :param X_b_i:
            :param y_i:
            :return:
            """
            return np.sum((X_b_k * (X_b_k.dot(theta) - y_k) ))* 2/len(X_b_k).

        def sgd(X_b, y, initial_theta, n_iters, t0=5, t1=50):

            def learning_rate(t):
                """
                计算学习率，t1 为了减慢变化速度，t0为了增加随机性
                :param t: 第t次循环
                :return:
                """
                return t0 / (t + t1)

            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                # 每次看k个元素
                i =0
                while i < m:
                    X_b_new = X_b[i:i+k]
                    y_new = y[i:i+k]
                    gradient = dJ_sgd(theta, X_b_new, y_new)
                    theta = theta - learning_rate(cur_iter * m + i+k) * gradient
                    i = i+k

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.random.randn(X_b.shape[1])
        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

7.2 随机

随机梯度下降法的优点

7.2-1

7.3 梯度上升法

7.3-1 7.3-2 7.3-3