线性回归原理、手写反向传播和pytorch实现

线性回归原理和实践

线性回归的数学推导

假设：我们知道机器学习的思想就是通过优化一个目标函数使其值达到最小，也就是找到函数的最小值点。线性回归的损失函数如下：
$$J(w)=\sum _{i=0}^n\frac{1}{2}(y_i-(w{x}_i+b){)}^2$$
这里的$w$和$b$是待学习得到的参数。下面从数学的角度推导一下为什么线性回归的目标函数会采用这样的公式？？
假设我们拿到一个数据集：

y	特征1	特征2	特征3
通过线性回归的方式拟合一系列的数据，数学公式如下：
$$y=w_1{x}_1+w_2{x}_2+w_3{x}_3+b$$
其中$w_1,w_2,w_3$分别表示不同特征对y值得影响程度，也就是待学习得参数，b在数学上表示偏置，也可以当成一个特征，只不过这列特征的数值全部是1.
如果我们的项目特征数量特别多，写成上面那样过于难看，所以数学上我们习惯写成矩阵相乘的形式，如下：
$$y=WX+b$$
其中W,X都是矩阵。
再假设我们用$\theta$来表示真实值与预测值之间的误差，
$$y_i=w_i{x}_i+b_i+{\theta }_i$$
这里的$i$代表每一组数据，$w_i,x_i$也是矩阵的形式（表示有多个特征）；
这里想必大家学过数学的都能看懂，那么下面就要用到大学的数学知识了。
假设1：所有的$\theta$都是服从独立同分布的；什么是独立，因为一个数据集的每条数据之间是相互独立的，没有任何关系。但这里你可能会想到时序数据，上一条数据会对下一条数据有一定的影响，那么还可以用这种方法吗，答案是肯定的，机器学习的思想就是基于独立同分布的，所以我们可以将时序数据中的前几个时刻的数据当成下一条数据的特征（这里和LSTM有点像，不再赘述了）。
假设2：这里的同分布是高斯分布（正态分布），记着是误差符合高斯分布，通过下图可以看出一般情况下落在$2\mu$之间，想一下误差呢也不会太大或者太小，几乎也都是落在这个区间上。
所以我们可以将上式写成：
高斯分布：$p({\theta }_i)=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{{\theta }_i^2}{2{\sigma }^2})$
代入：$p(y_i	x_i;w;b)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_i-(wx_i+b))2}{2\sigma2})$
这里引入似然函数，可能读者对似然函数有点忘了，这里我们简单的回忆一下，似然函数的定义：$L(x	\theta)=P(X=x	\theta)$,表示为在给定$x=\theta$的情况下，x变为$X$的概率。别忘了，我们这个高斯分布的含义，给定$x_i;w;b$的情况下，等于$y_i$（真实值）的概率，我们希望这个概率越大越好，灵光乍现，最大似然估计！
这里为了方便书写，将$w;b$合并为$w$。似然函数如下：
$$L(w)=\prod_{i=0}^np(y_i	x_i;w)=\prod_{i=0}n\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_i-(wTx_i))2}{2\sigma2})$$
为了方便求解，改为对数似然：
$$logL(w)=log\prod _{i=0}^n\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y_i-(w^T{x}_i){)}^2}{2{\sigma }^2})$$

$$logL(w)=\sum _{i=0}^{n}log\frac{1}{\sqrt{2\pi }\sigma }-\sum _{i=0}^n\frac{(y_i-(w^T{x}_i){)}^2}{2{\sigma }^2}$$
因为我们希望$logL(w)$越大越好，${\sum }_{i=0}^{n}log\frac{1}{\sqrt{2\pi }\sigma }$为常数，$\frac{1}{2{\sigma }^2}$也为常数，所以最后就希望${\sum }_{i=0}^n\frac{(y_i-(w^T{x}_i){)}^2}{2}$越小越好，也就是我们前面提到的损失函数。

下面利用偏导求解这个目标函数：

可以看出最后通过一系列的计算可以求出对应的权重，但是一般情况机器学习是求不出这个权重参数的，线性回归算是一个特例。

实践

下面通过实际的栗子来查看通过迭代梯度下降的方式和上面数学计算的结果进行比较。

import numpy as np
import matplotlib.pyplot as plt
def get_data(batch_size=8):
    # 产生随机数据，y = 3x+2
    np.random.seed(1)
    x = np.random.uniform(0,1,batch_size)*30
    y = 3*x + (1+np.random.random(batch_size))*2
    return x,y
x,y = get_data()
plt.scatter(x,y)
plt.show()

# 初始化权重
w = np.random.uniform(0,1,1)
b = np.random.uniform(0,1,1)
# print(w,b)
lr = 0.001
for i in range(10000):
    x,y = get_data()
    y_pred = w*x+b
    loss = 0.5*(y_pred - y)**2
    loss = sum(loss)
#     dw = x.dot((y_pred - y).T)
    dw = (y_pred-y).dot(x.T)
#     print(dw)
    db = sum(y_pred - y)
    w -=lr*dw
    b -=lr*db
    if i%100==0:
        x = np.arange(0,30)
        y = x*w + b
        plt.plot(x,y)
        x1,y1 = get_data()
        plt.scatter(x1,y1)
        plt.show()
        print(w,b)

w,b=[3.00987428] [2.87324644]

我们用上面求得的权重计算公式：$w=(X^{T}X{)}^{-}1X^{T}y$

x,y = get_data()
# 注意这里将x增加一组全1行，目的是求参数b
x = np.array([x, np.ones(len(x)])
res = np.matrix(x.dot(x.T)).I.dot(x).dot(y.T)
res
out:matrix([[3.00987428, 2.87324644]])

可以看出这种方法和迭代求梯度所求的结果是一样的，这也就验证了我们推导的过程是正确的。

用pytorch实现

import torch as t
from torch.autograd import Variable 
x,y = get_data()
w,b = np.random.uniform(0,1,1),np.random.uniform(0,1,1)
#将numpy转为tensor
x,y = Variable(t.Tensor(x)),Variable(t.Tensor(y))
w,b = Variable(t.Tensor(w),requires_grad=True),Variable(t.Tensor(b),requires_grad=True)
print(x,y,w,b)
out:tensor([1.2511e+01, 2.1610e+01, 3.4312e-03, 9.0700e+00, 4.4027e+00, 2.7702e+00,
        5.5878e+00, 1.0367e+01]) tensor([40.3255, 67.9068,  2.8487, 30.5804, 15.6169, 12.0667, 18.8182, 34.4414]) tensor([0.4173], requires_grad=True) tensor([0.5587], requires_grad=True)

lr=0.001
for i in range(1000):
#     x,y = get_data()
    y_pred = w.unsqueeze(0).mm(x.unsqueeze(0))+b
#     print(y_pred)
    loss = 0.5*(y_pred - y.unsqueeze(0))**2
#     print(loss)
    loss = loss.sum()
    loss.backward()
#     dw = x.dot((y_pred - y).T)
#     dw = (y_pred-y).mm(x.unsqueeze(0).t())
#     print(dw)
#     db = sum(y_pred - y)
    
    w.data.sub_(lr*w.grad.data)
    b.data.sub_(lr*b.grad.data)
    #梯度清零
    w.grad.data.zero_()
    b.grad.data.zero_()
   
    if i%100==0:
        x1 = np.arange(0,30)
        y1 = x1*w.data.numpy()[0] + b.data.numpy()[0]
        plt.plot(x1,y1)
        x2,y2 = get_data()
        plt.scatter(x2,y2)
        plt.show()
        print(w.data.numpy()[0],b.data.numpy()[0])

3.0101655 2.8694313

3.0100913 2.8704033

3.0100358 2.8711276

3.0099947 2.871668

3.009964 2.8720703

3.009941 2.8723693

3.0099242 2.8725924

3.0099115 2.8727586

3.0099018 2.8728845

3.0098948 2.8729756

以上就是线性回归的内容，后面还会进行更新。。。比如逻辑回归，决策树，随机森林家族等，支持向量机，等等