线性回归

最小二乘法
梯度下降法

最小二乘法

数学原理
- model: $y = \Sigma_{{i=1}}^{{n}}\omega_i x_i +b\\\ \ = \Sigma_{i = 0}^{n}\omega_ix_i (\omega_0=b,x_0=1)$
  
  m个sample，每个sample有n个特征，引入以下几个记号
  
  $X = \begin{bmatrix}x_{11}\ x_{12} \ ...x_{1n}\\x_{21}\ x_{22} \ ...x_{2n}\\...\ ...\ ...\ ...\\x_{m1}\ x_{m_2}...\ x_{mn}\end{bmatrix}=\begin{bmatrix}x_1^T\\x_2^T\\.\\.\\.\\x_m^T\end{bmatrix}\\ \omega = \begin{bmatrix}\omega_1\\\omega_2\\.\\.\\.\\\omega_n\end{bmatrix}\ y=\begin{bmatrix}y_1\\y_2\\.\\.\\.\\y_m\end{bmatrix}$
  
  我们有：
  
  $y = Xw$
- 均方差：MES $E(\omega)= \frac{1}{2m}\Sigma_{i = 1}^m(\hat y_i -y )^2$
  
  梯度： $\nabla E(\omega)=\frac{1}{m}\Sigma_{i=1}^m(\hat y_i -y )x_i\\\qquad =\frac{1}{m}X^T(X\omega-y)$
  - 注意：E的梯度是一个向量函数，是一个列向量
- 最小二乘：
  
  $\nabla E=0\\\Rightarrow X^T(X\omega-y)=0\\\Rightarrow\omega = (X^TX)^{-1}(X^Ty)$
在编程中遇到的问题
- python 面向对象编程-类的创建
  - 用class 创立类，调用类的方法之前先创建类的实例；默认类的构建函数init函数参数为空，如果自己定义init注意传参数
  - py和cpp不同，self相当于this指针；没有private，public声明，通过命名法区别类的私有函数和公有函数，代表是private的函数
- python导入数据集
  - 这里使用的是numpy.genfromtxt()函数
    
    genfromtxt(‘path’；delimiter’分割字符’还有很多参数)
- 用numpy处理数组
  - 矩阵切片：[ : , : ]“，”前表示行的切片
- 划分数据集
  - 用到sklearn.model_selection.train_test_split

import numpy as np

#最小二乘法
class olsLinearRegression:
    
    # X预处理
    def _preprocessing(self,X):
        m,n = X.shape
        X_ = np.empty((m,n+1))
        X_[:,1:] = X
        X_[:,0] = 1
        
        return X_
    
    #求self.w
    def _ols(self,X,y):
        temp = np.linalg.inv(np.matmul(X.T,X))
        temp_ = np.matmul(X.T,y)
        return np.matmul(temp,temp_)
    
    #训练model
    def train(self,X,y):
        X = self._preprocessing(X)
        self.w = self._ols(X,y)
    
    #做预测
    def predict(self,X):
        X = self._preprocessing(X)
        return np.matmul(X,self.w)

#导入数据集    
data = np.genfromtxt('E:\winequality-red.csv',delimiter = ';',skip_header = True)
X = data[:,:-1]
y = data[:,-1]

#划分数据集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size = 0.3)

#创建类的实例
ols = olsLinearRegression()
#训练模型并得到预测结果
ols.train(X_train,y_train)
y_pred = ols.predict(X_test)

#性能指标：均方差误差
from sklearn.metrics import mean_squared_error
mes = mean_squared_error(y_pred,y_test)
print(mes)

梯度下降

梯度下降的假设：凸函数，可微函数

不是凸函数存在局部最优解
不可微，无法获得封闭解

为什么有梯度下降

不是所有的问题都是凸优化问题，换而言之，不是所有的问题只有一个（极大/极小）值
当 $X^TX$ 不可逆， $\nabla E(\omega)$ 无零解，这时就要用到梯度下降了

梯度下降的概念

沿着某个点 $(x_i,y_i)$ 处梯度的反方向梯度下降
学习率：步长，
- 步长太小：收敛速度慢
- 步长太大：之字形收敛
$\omega := \omega-\eta \nabla E(\omega,x_i)$

隐含的问题

局部最优解
收敛速度
- $\nabla E 在各个\omega_i的分量大小不同，但乘的是相同的\eta,也就是说(\hat y_i-y)x_i越大，收敛速度越快$
- 因此可以通过数据的标准化 $\frac{x_{imax}-x_{ik}}{x_{imax}-x_{imin}}$ 减小这种影响，其中： $x_{imax}表示，样本x_i的各个特征数值最大的项；x_{ik}表示x_i的第k个特征的值$
- 也可以用变化学习率，选择最优学习率，这里不赘述，只呈现最基本的
  
  其实变化学习率的想法是很自然的，因为函数的梯度在不同点通常是不同的，梯度大的时候，相应的学习率应该小
$\omega$ 的初始化与样本点 $x_i$ 的选取
- $\omega的初始化这里不详细写$
- $x_i的选取，我们有随机梯度下降和批量梯度下降，小批量梯度下降$

批量梯度下降：

首先我们必须清楚梯度是怎么下降的

$\frac{\partial E}{\partial \omega_i}=0$

如果我们用的公式 $\omega := \omega-\eta \nabla E(\omega,x_i)\\\nabla E = \frac{1}{m}X^T(X\omega-y)$

梯度依赖于整个样本空间，那么我们称这种方法为：批量梯度下降法
一定能够求到求到全局最优解，但当样本数很多的时候，这种方法的效率比较低

随机梯度下降：

如果我们每步只使用一个样本来计算梯度；参数更新公式： $\omega := \omega - \eta(\omega^Tx_i-y_i)x_i$ ，梯度依赖于单个样本处的梯度
迭代一轮得到一个更新的参数向量，每轮迭代都要用到每一个样本点，适用于迭代次数远小于样本数的情况。注意不是每一轮迭代都是朝最优解方向走的，但总体趋势是趋近于最优解
原则上：校验每一轮更新结束后的 $\omega$ ，如果损失函数值稳定，训练结束，故需要设定阈值。又由于整体趋势是趋近于全局最优解的

一般是设置最大迭代次数（超参数）
思想：损失一部分精度（不适用精确的梯度），增加迭代次数；
评价：在非凸优化问题上，更容易跳过局部最优解，求得全局最优解；凸优化往往会损失一部分精度，但速度快

**小批量梯度下降 **

如果我们一次用k个样本点梯度的均值作为梯度来更新，那么这种方式就是小批量梯度下降
评价：介于批量梯度下降与随机梯度下降之间，同时具备二者的的优缺点，比较常用

 #基于梯度下降实现线性回归
 
import numpy as np

class GDLinearRegression:
    #创建构造函数
    def __init__(self,n_iter = 300,eta = 1e-2,tol = None):
        #训练迭代次数
        self.n_iter = n_iter
        #学习率
        self.eta = eta
        #误差变化阈值
        self.tol = tol
        #模型参数w(训练时初始化)
        self.w = None
        
    #预处理
    def _preprocessing(self,X):
       
        m,n = X.shape
        X_ = np.empty((m,n+1))
        X_[:,0] = 1
        X_[:,1:] = X
        
        return X_
    
    #损失函数
    def _loss(self,y,y_pred):
        
        return np.sum((y-y_pred)**2)/y.size
    
    #计算梯度
    def _gradient(self,X,y,y_pred):
        
        return np.matmul(y_pred-y,X)/y.size
    
    #梯度下降
    def _gradient_decent(self,w,X,y):
        
        #如果用户指定阈值，启用早期停止法
        if self.tol is not None:
            loss_old = np.inf
        #迭代次数
        for step_i in range(self.n_iter):
            #做预测
            y_pred =  self._predict(X,w)
            
            #计算损失
            loss = self._loss(y,y_pred)
            
            print("%4i Loss:%s"%(step_i,loss))
            
            #早期停止法
            if self.tol is None:
            # 如果损失值小于阈值终止迭代
                if loss_old - loss <self.tol:
                    break
                less_old = loss
            
            #计算梯度
            grad = self._gradient(X,y,y_pred)
            w -= self.eta*grad
            
        
    #做训练
    def train(self,X_train,y_train):
        X_train = self._preprocessing(X_train)
        #初始化参数向量w//这里_作为一个临时的名字，后面不会在用这个变量
        _,n = X_train.shape
        self.w = np.random.random(n)*0.05
        #执行梯度下降训练
        self._gradient_decent(self.w,X_train,y_train)
        
    #做预测的内部接口函数//对data的每个数据都进行预测
    def _predict(self,X,w):
        
        return np.matmul(X,w)
    
    #做预测
    def predict(self,X):
        X = self._preprocessing(X)
        return self._predict(X,self.w)
    

#导入数据集
data = np.genfromtxt('E:\winequality-red.csv',delimiter = ';',skip_header = True)
X = data[:,:-1]
y = data[:,-1]

#划分数据集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size = 0.3)

#对X进行列标准化，优化梯度下降
from sklearn.preprocessing import StandardScaler
ss = StandardScaler()

'''
#计算特征的均值和标准差
ss.fit(X_train)
StandardScaler(copy = True,with_mean = True,with_std = True)

'''

#标准化
X_train_std = ss.fit_transform(X_train)
X_test_std = ss.fit_transform(X_test)


#创建类的实例
GDS = GDLinearRegression(n_iter = 3700,eta = 0.01,tol = 0.000001)

#训练模型
GDS.train(X_train_std,y_train)
y_pred = GDS.predict(X_test_std)

from sklearn.metrics import mean_squared_error
mes = mean_squared_error(y_pred,y_test)
ptint(mes)

将两个模型放在一起，更容易看到梯度下降的优化效果

import numpy as np

#最小二乘法
class olsLinearRegression:
    
    # X预处理
    def _preprocessing(self,X):
        m,n = X.shape
        X_ = np.empty((m,n+1))
        X_[:,1:] = X
        X_[:,0] = 1
        
        return X_
    
    #求self.w
    def _ols(self,X,y):
        temp = np.linalg.inv(np.matmul(X.T,X))
        temp_ = np.matmul(X.T,y)
        return np.matmul(temp,temp_)
        
    def train(self,X,y):
        X = self._preprocessing(X)
        self.w = self._ols(X,y)
    
    def predict(self,X):
        X = self._preprocessing(X)
        return np.matmul(X,self.w)

class GDLinearRegression:
    #创建构造函数
    def __init__(self,n_iter = 200,eta = 1e-2,tol = None):
        #训练迭代次数
        self.n_iter = n_iter
        #学习率
        self.eta = eta
        #误差变化阈值
        self.tol = tol
        #模型参数w(训练时初始化)
        self.w = None
        
    #预处理
    def _preprocessing(self,X):
       
        m,n = X.shape
        X_ = np.empty((m,n+1))
        X_[:,0] = 1
        X_[:,1:] = X
        
        return X_
    
    #损失函数
    def _loss(self,y,y_pred):
        
        return np.sum((y-y_pred)**2)/y.size
    
    #计算梯度
    def _gradient(self,X,y,y_pred):
        
        return np.matmul(y_pred-y,X)/y.size
    
    #梯度下降
    def _gradient_decent(self,w,X,y):
        
        #如果用户指定阈值，启用早期停止法
        if self.tol is not None:
            loss_old = np.inf
        #迭代次数
        for step_i in range(self.n_iter):
            #做预测
            y_pred =  self._predict(X,w)
            
            #计算损失
            loss = self._loss(y,y_pred)
            
            print("%4i Loss:%s"%(step_i,loss))
            
            #早期停止法
            if self.tol is None:
            # 如果损失值小于阈值终止迭代
                if loss_old - loss <self.tol:
                    break
                less_old = loss
            
            #计算梯度
            grad = self._gradient(X,y,y_pred)
            w -= self.eta*grad
            
        
    #做训练
    def train(self,X_train,y_train):
        X_train = self._preprocessing(X_train)
        #初始化参数向量w//这里_作为一个临时的名字，后面不会在用这个变量
        _,n = X_train.shape
        self.w = np.random.random(n)*0.05
        #执行梯度下降训练
        self._gradient_decent(self.w,X_train,y_train)
        
    #做预测的内部接口函数//对data的每个数据都进行预测
    def _predict(self,X,w):
        
        return np.matmul(X,w)
    
    #做预测
    def predict(self,X):
        X = self._preprocessing(X)
        return self._predict(X,self.w)
    

data = np.genfromtxt('E:\winequality-red.csv',delimiter = ';',skip_header = True)
X = data[:,:-1]
y = data[:,-1]

#划分数据集
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size = 0.3)

#创建类的实例
ols = olsLinearRegression()
ols.train(X_train,y_train)
y_pred = ols.predict(X_test)


#对X进行列标准化，优化梯度下降
from sklearn.preprocessing import StandardScaler
ss = StandardScaler()

'''
#计算特征的均值和标准差
ss.fit(X_train)
StandardScaler(copy = True,with_mean = True,with_std = True)

'''

#标准化
X_train_std = ss.fit_transform(X_train)
X_test_std = ss.fit_transform(X_test)


#创建类的实例
GDS = GDLinearRegression(n_iter = 3700,eta = 0.001,tol = 0.000001)

#训练模型
GDS.train(X_train_std,y_train)
y_pred_gds = GDS.predict(X_test_std)

#性能指标：均方差误差
from sklearn.metrics import mean_squared_error
mes = mean_squared_error(y_pred,y_test)
print("OLS mes ")
print(mes)
mes_ = mean_squared_error(y_pred_gds,y_test)
print(mes_)

LinearRegression_补充了梯度下降的完整版

LinearRegression_补充了梯度下降的完整版

线性回归

最小二乘法

梯度下降

将两个模型放在一起，更容易看到梯度下降的优化效果

相关阅读更多精彩内容

友情链接更多精彩内容