机器学习基础：线性回归小结

机器学习基础：线性回归小结

CONTENT

1.线性回归形式和损失函数

• 基本形式
• 损失函数

2.线性回归的求解算法

• 梯度下降法
• 最小二乘法
• 牛顿法
• 拟牛顿法

3.多项式回归和广义线性回归

• 多项式回归
• 广义线性回归

4.线性回归中的正则化

• L1正则： Lasso
• L2正则： Ridge

5.代码实践

• 数据准备
• sklearn 线性回归
• scripts with dirty hands
  • 梯度下降法求解
  • 最小二乘法求解
  • 牛顿法求解
  • Ridge回归

1. 线性回归形式和损失函数

1.1 基本形式

\[h_θ(x_1,x_2,...x_n)=θ_0+θ_1x_1+...+θ_nx_n \]

简化表示为

\[h_θ(x_0,x_1,...x_n)=\sum_{i=0}^{n}θ_ix_i \]

1.2 损失函数

代数法表示：

\[J(θ_0,θ_1...,θ_n)=\sum_{j=1}^{m}(\sum_{i=0}^{n}\theta_ix_i-y_j) \]

矩阵表示：

\[J(θ)=\frac{1}{2}(Xθ−Y)^T(Xθ−Y) \]

注： m为样本个数， n为特征个数加一

2. 线性回归的求解算法

2.1 梯度下降算法

设定初始参数,不断迭代，使得最小化，θ的迭代公式如下：

\[θ=θ−αX^T(Xθ−Y) \]

当J为凸函数时，梯度下降法相当于让参数\(\theta\)不断向J的最小值位置移动
梯度下降法的缺陷：如果函数为非凸函数，有可能找到的并非全局最优值，而是局部最优值。

2.2 最小二乘法

令：

\[\frac{\partial}{\partial\theta}J(\theta) = X^T(X\theta-Y) = 0 \]

利用矩阵的链式求导法则，和两个矩阵的求导公式：

\[公式一： \partial_{X}(X^TX) = 2X \\ 公式二： \nabla_{X}f(AX+B) = A^T\nabla_{Y}f, Y = AX+B \]

整理得到：

\[θ=(X^TX)^{−1}X^TY \]

1. 最小二乘法需要计算\(X^TX\)的逆矩阵，有可能它的逆矩阵不存在，这样就没有办法直接用最小二乘法了，此时梯度下降法仍然可以使用。当然，我们可以通过对样本数据进行整理，去掉冗余特征。让\(X^TX\)的行列式不为0，然后继续使用最小二乘法。
2. 当样本特征n非常的大的时候，计算\(X^TX\)的逆矩阵是一个非常耗时的工作（nxn的矩阵求逆），甚至不可行。此时以梯度下降为代表的迭代法仍然可以使用。那这个n到底多大就不适合最小二乘法呢？如果你没有很多的分布式大数据计算资源，建议超过10000个特征就用迭代法吧。或者通过主成分分析降低特征的维度后再用最小二乘法。
3. 如果拟合函数不是线性的，这时无法使用最小二乘法，需要通过一些技巧转化为线性才能使用，此时梯度下降仍然可以用
4. 特殊情况下，当样本量m很少，小于特征数n的时候，这时拟合方程是欠定的，常用的优化方法都无法去拟合数据。当样本量m等于特征数n的时候，用方程组求解就可以了。当m大于n时，拟合方程是超定的，也就是我们常用与最小二乘法的场景了。

2.3 牛顿法

将\(J(\theta)\)泰勒展开到二阶：

\[J(\theta) = J(\theta_0) + J^{'}(\theta_0)(\theta-\theta_0) + \frac{1}{2}J^{''}(\theta_0)(\theta-\theta_0)^2 \]

两边求导，得：

\[\theta = \theta_0 - \frac{J^{'}(\theta_0)}{J^{''}(\theta_0)} \]

重复迭代式可以利用Hessian矩阵：

\[\theta =\theta - H^{-1}\Delta_{\theta}l(\theta) \]

牛顿法的收敛速度非常快，但海森矩阵的计算较为复杂，尤其当参数的维度很多时，会耗费大量计算成本。可以用其他矩阵替代海森矩阵，用拟牛顿法进行估计

2.4 拟牛顿法

拟牛顿法的思路是用一个矩阵替代计算复杂的海森矩阵H，因此要找到符合H性质的矩阵。

要求得海森矩阵符合的条件，同样对泰勒公式求导\(f'(x) = f'(x_0) + f''(x_0)x -f''(x_0)x_0\)
令\(x = x_1\)，即迭代后的值，代入可得：

\[f'(x_1) = f'(x_0) + f''(x_0)x_1 - f''(x_0)x_0 \]

更一般的

\[𝑓′(𝑥𝑘+1)=𝑓′(𝑥𝑘)+𝑓″(𝑥𝑘)𝑥𝑘+1−𝑓″(𝑥𝑘)𝑥𝑘f′(xk+1)=f′(xk)+f″(xk)xk+1−f″(xk)xk \\ 𝑓′(𝑥𝑘+1)−𝑓′(𝑥𝑘)=𝑓″(𝑥𝑘)(𝑥𝑘+1−𝑥𝑘)=𝐻(𝑥𝑘+1−𝑥𝑘)f′(xk+1)−f′(xk)=f″(xk)(xk+1−xk)=H(xk+1−xk) \]

x_k为第k个迭代值

即找到矩阵G，使得它符合上式。 常用的拟牛顿法的算法包括DFP，BFGS等。

3. 多项式回归和广义线性回归

3.1 多项式回归

回到开始的线性模型

\[h_\theta(x_1,x_2,...x_n)=θ_0+θ_1x_1+...+θ_nx_n \]

, 如果这里不仅仅是x的一次方，比如增加二次方，那么模型就变成了多项式回归。这里写一个只有两个特征的p次方多项式回归的模型：

\[h_θ(x_1,x_2)=θ_0+θ_1x_1+θ_2x_2+θ_3x^2_1+θ_4x^2_2+θ_5x_1x_2 \]

令\(*x*_0=1,x_1=x_1,x_2=x_2,x_3=x^2_1,x_4=x^2_2,x_5=x_1x_2\)

这样我们就得到了下式：

\[h_θ(x_1,x_2)=θ_0+θ_1x_1+θ_2x_2+θ_3x_3+θ_4x_4+θ_5x_5 \]

可以发现，又重新回到了线性回归，这是一个五元线性回归，可以用线性回归的方法来完成算法。对于每个二元样本特征(x1,x2), 得到一个五元样本特征(1,x1,x2,x21,x22,x1x2)，通过这个改进的五元样本特征，重新把不是线性回归的函数变回线性回归。

3.2 广义线性回归

输出Y不满足和X的线性关系，但是lnY 和X满足线性关系，模型函数如下：

\[\ln{Y} = X\theta \]

这样对与每个样本的输入y，我们用 lny去对应， 从而仍然可以用线性回归的算法去处理这个问题。

4. 线性回归中的正则化

4.1 L1正则

Lasso 回归

\[J(θ)=\frac{1}{2}(Xθ−Y)^T(Xθ−Y) + \alpha\left\|\theta\right\|_1 \]

Lasso回归可以使得一些特征的系数变小，甚至还是一些绝对值较小的系数直接变为0。增强模型的泛化能力。Lasso回归的求解办法一般有坐标轴下降法（coordinate descent）和最小角回归法（ Least Angle Regression）。

4.2 L2正则

Ridge回归

\[J(θ)=\frac{1}{2}(Xθ−Y)^T(Xθ−Y) + \frac{1}{2}\alpha\left\|\theta\right\|_2^2 \]

Ridge回归在不抛弃任何一个特征的情况下，缩小了回归系数，使得模型相对而言比较的稳定，但和Lasso回归比，这会使得模型的特征留的特别多，模型解释性差。

Ridge回归的求解比较简单，一般用最小二乘法。这里给出用最小二乘法的矩阵推导形式，和普通线性回归类似。

\[\theta = (X^TX + \alpha E)^{-1}X^TY \]

5. 代码实践

5.1 数据准备

view code#生成数据
import numpy as np
#生成随机数
np.random.seed(1234)
x = np.random.rand(500,3)
#构建映射关系，模拟真实的数据待预测值,映射关系为y = 4.2 + 5.7*x1 + 10.8*x2
y = x.dot(np.array([4.2,5.7,10.8]))

5.2 sklearn 线性回归

view codeimport numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
%matplotlib inline

# 调用模型
lr = LinearRegression(fit_intercept=True)
# 训练模型
lr.fit(x,y)
print("估计的参数值为：%s" %(lr.coef_))
# 计算R平方(拟合优度，取值[0,1], 越接近1拟合越好)
print('R2:%s' %(lr.score(x,y)))
# 任意设定变量，预测目标值
x_test = np.array([2,4,5]).reshape(1,-1)
y_hat = lr.predict(x_test)
print("预测值为: %s" %(y_hat))

5.3 scripts with dirty hands

5.3.1 梯度下降法求解

view codeclass LR_GD():
    def __init__(self):
        self.w = None     
    def fit(self,X,y,alpha=0.02,loss = 1e-10): # 设定步长为0.002,判断是否收敛的条件为1e-10
        y = y.reshape(-1,1) #重塑y值的维度以便矩阵运算
        [m,d] = np.shape(X) #自变量的维度
        self.w = np.zeros((d)) #将参数的初始值定为0
        tol = 1e5
        while tol > loss:
            h_f = X.dot(self.w).reshape(-1,1) 
            theta = self.w + alpha*np.mean(X*(y - h_f),axis=0) #计算迭代的参数值
            tol = np.sum(np.abs(theta - self.w))
            self.w = theta
    def predict(self, X):
        # 用已经拟合的参数值预测新自变量
        y_pred = X.dot(self.w)
        return y_pred  

if __name__ == "__main__":
    lr_gd = LR_GD()
    lr_gd.fit(x,y)
    print("估计的参数值为：%s" %(lr_gd.w))
    x_test = np.array([2,4,5]).reshape(1,-1)
    print("预测值为：%s" %(lr_gd.predict(x_test)))

5.3.2 最小二乘法求解

view codeclass LR_LS:
    def __init__(self, fit_intercept=True):
        self.w = None
        self.fit_intercept = fit_intercept # bias 截距

    def fit(self, X, y):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]

        pseudo_inverse = np.dot(np.linalg.inv(np.dot(X.T, X)), X.T)
        self.w = np.dot(pseudo_inverse, y)

    def predict(self, X):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]
        return np.dot(X, self.w)

5.3.3 牛顿法求解

步骤如下：

view codedef first_derivativ(X, y, w):
    return np.matmul(X.T, np.matmul(X, w) - y)

def second_derivative(X):
    return np.matmul(X.T, X)

def get_error(X, y, w):
    tmp = np.matmul(X, w) - y
    return np.matmul(tmp.T, tmp) / 2
 
def get_min_m(X, y, sigma, delta, d, w, g):
    m = 0
    while True:
        w_new = w + pow(sigma, m) * d
        left = get_error(X, y , w_new)
        right = get_error(X, y , w) + delta * pow(sigma, m) * g.T * d
        if left <= right:
            break
        else:
            m += 1
    return m   

class LR_NEWTON():
    def __init__(self):
        self.w = None
    
    def fit(self, X, y, max_iter=50, sigma=0.5, delta=0.25):
        '''
            sigma : (0, 1)
            delta : (0. 0.5)
        '''
        print(type(X))
        n = X.shape[1]
        self.w = np.mat(np.ones((n, 1)))
        it = 0
        while it <= max_iter:
            g = first_derivativ(X, y, self.w)
            G = second_derivative(X)
            try:
                d = -np.linalg.inv(G) * g
            except Exception:
                print('不可逆')
                break
            m = get_min_m(X, y, sigma, delta, d, self.w, g)  # 得到最小的m
            self.w = self.w + pow(sigma, m) * d
            if it % 10 == 9:
                print ("\t---- itration: ", it, " , error: ", get_error(X, y , self.w)[0, 0])
            it += 1
    def predict(self, X):
        y_pred = X.dot(self.w)
        return y_pred 
    
if __name__ == "__main__":
    lr_nt = LR_NEWTON()
    lr_nt.fit(x,y.reshape(-1,1))
    print("估计的参数值为：%s" %(lr_nt.w))
    x_test = np.array([2,4,5,2,4,5,2,4,5]).reshape(1,-1)
    print("预测值为：%s" %(lr_nt.predict(x_test)))

5.3.4 Ridge回归

view codeclass RidgeRegression:
    def __init__(self, fit_intercept=True):
        self.w = None
        self.fit_intercept = fit_intercept

    def fit(self, X, y, alpha=1):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]

        A = alpha * np.eye(X.shape[1])
#         A = 0
        pseudo_inverse = np.dot(np.linalg.inv(np.dot(X.T, X) + A), X.T)
        self.w = np.dot(pseudo_inverse, y)

    def predict(self, X):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]
        return np.dot(X, self.w)
    
if __name__ == "__main__":
    lr_ridge = RidgeRegression()
    lr_ridge.fit(x,y.reshape(-1,1))
    print("估计的参数值为：%s" %(lr_ridge.w))
    x_test = np.array([2,4,5,2,4,5,2,4,5]).reshape(1,-1)
    print("预测值为：%s" %(lr_ridge.predict(x_test)))