【机器学习】线性回归理论和实践

线性回归代码详解

1. 理论：

模型：线性模型；策略：均方误差最小化；算法：梯度下降算法

定义：给定数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x n , y n ) } D=\{(\bm{x}_1,y_1),(\bm{x}_2, y_2),\cdots,(\bm{x}_n,y_n)\} D={(x1​,y1​),(x2​,y2​),⋯,(xn​,yn​)}，其中 ${\mathbf{x}}_i=(x_{i1};x_{i2};\cdots ;x_{id})$，$d$ 是特征维度，$y\in \mathbb{R}$. 线性回归就是寻找函数 $f({\mathbf{x}}_i)=y_i$，即使得预测值等于输出值. 对于某一个样本${\mathbf{x}}_i$的模型输出表示.

模型
$$f({\mathbf{x}}_i)=b+\sum _i^d{w}_i{\mathbf{x}}_i\text{\hspace{1em}(1)}$$

策略
$$f^{\ast }=\arg \underset{f}{\min }L(f)\text{\hspace{1em}(2)}$$

$${\mathbf{w}}^{\ast },b^{\ast }=\arg \underset{\mathbf{w},b}{\min }L(\mathbf{w},b)=\arg \underset{\mathbf{w},b}{\min }\sum _{i=1}^n(y_i-(b+\mathbf{w}{\mathbf{x}}_i){)}^2\text{\hspace{1em}(3)}$$

2. 代码实现

2.1 梯度下降算法求解

算法（Gradient Descent）—前提L是可微分的

流程：参数 $\mathbf{w}$，学习率为 $\eta$

1. 随机初始化 ${\mathbf{w}}^0$;
2. 计算梯度 $\frac{dL}{d\mathbf{w}}{|}_{\mathbf{w}={\mathbf{w}}^0}$;
3. 更新梯度 ${\mathbf{w}}^1\leftarrow {\mathbf{w}}^0-\eta \frac{dL}{d\mathbf{w}}{|}_{\mathbf{w}={\mathbf{w}}^1}$;
4. 重复步骤2~3更新参数

每个参数都进行GD寻找 local optimal，对于 Linear Regression 可以得到最优解，因为 L 是 convex.

求解：

$$\frac{\partial L}{\partial w}=\frac{1}{N}\sum _{i=1}^{n}2(y_i-(b+w\cdot x_i))(-x_i)\text{\hspace{1em}(4)}$$

$$\frac{\partial L}{\partial b}=\frac{1}{N}\sum _{i=1}^{n}2(y_i-(b+w\cdot x_i))(-1)\text{\hspace{1em}(5)}$$

选择模型的时候，可以选择其他模型，比如引入 $x$ 高次项，即选择更加复杂的模型，可能会造成过拟合，因为模型比较复杂，训练集强行拟合使得测试集效果并不好.

多个线性组合：引入$\delta$，针对不同类型设置不同的线性模型。

正则化：$L={\sum }_n(y_i-(b+\sum w_i{x}_i){)}^2+\lambda \sum (w_i{)}^2$，加上正则项使得函数更加平滑，不考虑bias，bias对bias平滑程度没有影响

import numpy as np

def stepGradient(bias, weight, data, learning_rate, iterations):
    """
    Args:
        b_current: (1, 1)
        w_current: (1, d)
        data: (n, d)
        learning_rate: R
    """
    # 梯度清零
    w_grad = np.zeros(data.shape[1]-1)          # (m, 1)
    b_grad = 0.0
    # 得到样本数据
    bias = np.mat(bias)                         # (1, 1)
    weight = np.mat(weight)                     # (1, m)
    x = np.mat(data[:, 0:-1])                   # (n, m)
    y = np.mat(data[:, -1].reshape([-1, 1]))    # (n, 1)
    # 样本数
    N = float(len(data))

    for j in range(iterations):
        f_x = x * weight.T + bias     # [n, 1]
        
        loss = (y - f_x).T * (y - f_x) / N
        w_grad = - (2/N) * ((y - f_x).T * x)
        b_grad = - (2/N) * ((y - f_x).T * np.ones(shape=[x.shape[0],1]))
        # for i in range(len(data)):
        #     # 计算梯度之和
        #     w_grad = w_grad - (2.0/N) * (y[i]- bias - weight * x[i].T) * x[i]
        #     b_grad = b_grad - (2.0/N) * (y[i]- bias - weight * x[i].T)

        # 更新梯度
        weight = weight - learning_rate * w_grad
        bias = bias - learning_rate * b_grad
        # if j % 50 == 0:
            # print(loss)

    return bias, weight

# 数据准备
x_train = np.mat(np.random.normal(0,10,[5,3]))
Weights=np.mat(np.array([[3,4,6]]))
y_train = np.dot(x_train, Weights.T)+5
data = np.concatenate((x_train, y_train), axis=1)              # data: [n, m] 

# 参数
b = np.mat(0.0)                                                 # (1, 1)
w = np.mat(np.zeros(data.shape[1]-1))                           # (1, m)        # 使用行向量
lr = 0.001
iterations = 10000

b, w = stepGradient(b, w, data, lr, iterations)                 
pred = b + x_train * w.T                                       # 预测结果

print(w)
print(y_train)
print(pred)

线性回归需清楚推导的损失对参数的偏导数，得到梯度后进行参数的更新，注意矩阵乘法中的细节，w的形状是什么，b的形状是什么，在整个算法中，一行为一个样本，一列为一个特征随机变量。清楚了这些后便可以得到使用梯度下降算法进行求解

2.2 最小二乘法求解

参考西瓜书，数据集D: n x (m+1)矩阵 $\mathbf{X}$，输入输出对于单个样本都是列向量格式，优化函数是：

$${\hat{\mathbf{w}}}^{\ast }=\underset{\hat{w}}{\mathrm{argmin}}(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}{)}^{\mathrm{T}}(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})\text{\hspace{1em}(1)}$$

注意这里数据矩阵，n为样本数，m为样本特征维度，与西瓜书表达有所区别，为了方便学习，都将一行作为一个样本，一列作为一个特征维度。

令$E_{\hat{\mathbf{w}}}=(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}{)}^{\mathrm{T}}(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})$，对 $\hat{\mathbf{w}}$ 求导得到

$$\frac{\partial E_{\hat{\mathbf{w}}}}{\partial \hat{\mathbf{w}}}=2{\mathbf{X}}^{\mathrm{T}}(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y})\text{\hspace{1em}(2)}$$

当${\mathbf{X}}^{\mathrm{T}}\mathbf{X}$ 为满秩矩阵或正定矩阵时，令上式为零可得

$$\widehat{{\mathbf{w}}^{\ast }}=({\mathbf{X}}^{\mathrm{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}\text{\hspace{1em}(3)}$$

最后线性回归矩阵为

$$f(\widehat{{\mathbf{x}}_i})={\hat{\mathbf{x}}}_i^{\mathrm{T}}({\mathbf{X}}^{\mathrm{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}\text{\hspace{1em}(4)}$$

import numpy as np

# 准备数据
x_train = np.random.normal(0,10,[5,3])               # [n, m]
# print(x_train.shape)
Weights = np.array([[3,4,6]])                        # [1, m]
y_train = np.dot(x_train, Weights.T).squeeze()+5     # [1, n]

b = np.ones((x_train.shape[0], 1))                   # [n, 1]
X = np.concatenate([x_train, b], axis=1)             # [n, m]

# 检查矩阵是否满秩
rank = np.linalg.matrix_rank(X)
print(rank)

# 直接求解闭式解得到参数结果
w_ = (np.mat(X).T*np.mat(X)).I*(np.mat(X).T)*(np.mat(y_train).T)

# 预测输出
pred = X*w_
print(w_)
print(y_train)
print(pred)