实现梯度下降（线性回归为例）

2021-01-02

梯度下降（线性回归为例）

0. 相关工作

0.1 线性回归模型

$\hat {y} = \theta_0 + \theta_1x_1 + \theta_2x_2 +···+\theta_nx_n$

其中，$n$ 为 $x$ 的特征数量。

转为矩阵相乘形式：

$y = X·\theta$

其中，$X$ 为样本矩阵，个人喜欢行数代表样本数量 $m$，列数代表特征维度 $n$。$\theta$ 为参数矩阵，大小为 $n * 1$。

若行数代表特征维度 $n$，列数代表样本数量 $m$，可写为：

$y = \theta^{T}·X$

0.2 定义损失函数

采用均方误差损失函数$Mean Square Error(MSE) $ 。

某个样本的损失函数定义如下：

$loss_j = (X_{j}·\theta - y)^2$

整个训练集的损失函数为：

$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m} (X_i·\theta - y_i)^2$

我们要求的就是使得$J(\theta)$最小的$\theta$。

$\hat{\theta} = arg\min_{\theta} J(\theta)$

0.3 定义梯度公式

梯度就是由多维变量偏导数的向量。

可知$J(\theta)$ 是一个复合函数，求导时采用链式法则，$J$ 对每个维度的参数的偏导，定义为：

$\frac{\partial J}{\partial \theta_j} = \frac{1}{m} \sum_{i=1}^{m}(X^{i}·\theta - y^i)X_{j}^{i}$

其中$X_j^i$是第$j$个属性之前的系数。表示第i个样本第j个维度的特征值。

矩阵相乘形式：

$\frac{\partial J}{\partial \theta_j} = \frac{1}{m} X_j^T(X·\theta - y)$

$X_j^T$表示样本矩阵第j个维度的所有特征值。

由此进一步推导出：

$\frac{\partial J}{\partial \theta} = \frac{1}{m} X^T(X·\theta - y)$

1. 批量梯度下降

$\theta_{t+1} = \theta_t - \alpha\nabla\theta$

import numpy as np
from matplotlib import pyplot as plt
#定义损失函数
def loss(W,x,y):
    m = len(y)
    predictions = np.dot(x,W)
    return (1/2*m)*np.sum(np.square(predictions - y))

#定义批量梯度公式
def gradient(W,x,y):
    m = len(y)
    predictions = np.dot(x,W)
    return (1/m)*x.T.dot(predictions-y)
#迭代训练
def gradient_descent(W,x,y,alpha=0.001,iterations=1000):
    cost_history = np.zeros(iterations)
    
    for i in range(iterations):
        cost_history[i] = loss(W,x,y)
        grd = gradient(W,x,y)
        print(grd.shape)#(x.shape[1],1)
        W = W - alpha*grd
    return W,cost_history
if __name__ == "__main__":
    #生成数据
    X = 2 * np.random.rand(100,1)
    y = 10 +3 * X + np.random.randn(100,1)

    W = np.random.randn(2,1)
    #偏置项
    X_b = np.ones((X.shape[0],1))
    #加入偏置项，偏执项为(100,1)的全为1的向量
    X_ = np.hstack((X_b,X))#(100,2)
    iteration = 10000
    W,loss = gradient_descent(W,X_,y,alpha=0.01,iterations=iteration)

    print(W)
    #由于将偏置项放在X的第一列，也就是第0维。因此，bias = 8.67 weight = 3.85
    #array([[8.67047187],
    #      [3.85509216]])

    y_predict = np.dot(X_,W)
    it = np.linspace(1,iteration,iteration)
    plt.plot(it,loss)

    plt.scatter(X,y)
    plt.plot(X,y_predict)

2. 随机梯度下降

随机梯度下降是指，在进行梯度更新时，随机选取某一个样本来更新梯度，而非对整个样本数据集求梯度。

def stochastic_gradient_descent(W,x,y,alpha,iterations = 1000):
    m = len(y)
    cost_history = np.zeros(iterations)
    for i in range(iterations):
        predictions = np.dot(x,W)
        cost_history[i] = (1/2*m)*np.sum(np.square(predictions - y))
        
        rand_index = np.random.randint(0,m)
        #从数据集中取出索引为rand_index的数据
        x_i = x[rand_index,:].reshape(1,x.shape[1])
        y_i = y[rand_index].reshape(1,1)
        grad = gradient(W,x_i,y_i)
        
        W = W - alpha*grad
    return W,cost_history

3. 小批量梯度下降(mini-batch)

小批量梯度下降指：每一步的梯度计算，既不是基于整个训练集（如批量梯度下降）也不是基于单个实例（如随机梯度下降），而是基于一小部分随机的实例集也就是小批量，来进行梯度的更新。

def mini_batch_gradient_descent(W,x,y,alpha=0.01,itera=1000,batch_size=10):
	m = len(y)
    cost_history = []
    batches = int(m / batch_size)
    
    for it in range(itera):
       	cost = 0.0
        
        for i in range(0,m,batch_size):
            X_i = x[i:i+batch_size]#size:[batch_size,m]
            y_i = y[i:i+batch_size]#size:[batch_size,1]
        	
            predictions = np.dot(X_i,W)
            #计算整个数据集上的loss
            cost += (1/2*m)*np.sum(np.square(predictions - y_i))
            grad = gradient(W,X_i,y_i)
            
            W = W - alpha*grad
        cost_history.append(cost)
    return W,cost_history

4. 动量梯度下降

$v_t = \gamma v_{t-1} + \alpha \nabla J(\theta)$ $\theta = \theta - v_t$

$\alpha$ 即学习率。$\gamma$ 是动量系数。

def momentum_gradient_descent(W,x,y,alpha=0.01,mini_batch=20,itera=1000):
    volocity = np.zeros(W.shape)
    gamma = 0.9
    m = len(y)
    cost_history = []
    for it in range(itera):
        cost = 0.0
        
        for i in range(0,m,batch_size):
            X_i = x[i:i+batch_size]#size:[batch_size,m]
            y_i = y[i:i+batch_size]#size:[batch_size,1]

            predictions = np.dot(X_i,W)
            #计算整个数据集上的loss
            cost += (1/2*m)*np.sum(np.square(predictions - y_i))
            grad = gradient(W,X_i,y_i)
            #volocity的维度和W相同，列向量 直接矩阵计算
            #相当于在每个维度上计算gamma*volocity[dim] + alpha*grad[dim]
            volocity = gamma * volocity + alpha*grad 
            W = W - volocity
        cost_history.append(cost)
    return W,cost_history