线性回归_输入x有m个样本,每个样本有n个特征,写出一下变量的纬度-CSDN博客

本文链接：https://blog.csdn.net/gadwgdsk/article/details/79636147

1、心路历程

近来学习了感知器、线性回归、逻辑回归等，发现他们大体相似，线性回归其实是做数据拟合，而感知器可以认为是将线性回归应用到分类去了，使用了符号函数来离散化结果，而逻辑回归本质上也是寻找一维矩阵，只是它对线性方程用了sigmod函数，这个方法的好处是可以得到样本属于某一类的概率。

线性回归一般使用最小二乘法来实现，常见有两种推导形式，一种是适用于简单二维平面的场景，y=wx+b，推导过程比较简单，而对于多维场景则需要使用矩阵推导，它适用于所有场景。但是对于高维度、大数据集的回归求解，直接使用矩阵计算也不适用，要使用梯度下降法来做，代码实现中主要以梯度下降法为主。

2、一元线性回归

也就是

这里就是要求w和b

均方差是常用的性能度量方法，尝试另f(xi)与yi的均方差最小，于是有

（至于为什么最小二乘法是最优目标函数，可以使用极大似然估计来推导。由于误差是多个独立的因素共同影响得到的，根据中心极限定理，可以认为其服从高斯分布，参考：https://www.cnblogs.com/futurehau/p/6105011.html）

因此在线性回归中，最小二乘法（Least Mean Squre，LMS）就是试图找到一条直线，是的样本到直线的欧式距离之和最小，求解w和b就是目标函数最小化的过程，其过程称为线性回归的最小二乘参数估计（parameter estimation）。由于目标函数在凸函数，因此目标函数有最小值。判断函数是否凸函数方法如下：

目标函数的导数为：

其二阶导显示都是非负的。另二者导数为零联立求解可以得到：

由于上面的变量都是标量，因此计算很简单。

3、多元线性回归

这里吧w和b合并起来了，用w0表示b

我们需要最小化目标函数，关心theta取什么值的时候，目标函数取得最小值，而目标函数连续，那么theta一定为目标函数的驻点，所以我们求导寻找驻点

求导可得：

当是满秩矩阵（满秩矩阵是判断矩阵是否可逆的充要条件）或正定矩阵(若A是正定矩阵，必有|A|>0，则A必可逆)时候，令上式为0可以得到：

若特征数比样本数量还多，

不是满秩矩阵，无法求逆矩阵，不过这种情况现实中好像比较少见，一般样本数据都比较多。

关于矩阵的求导，维基百科写的很详细，可以查表

https://en.wikipedia.org/wiki/Matrix_calculus#Scalar-by-vector_identities

上述最后一步有一些问题，假如X'X不可逆呢？

我们知道X'X是一个办正定矩阵，所以若X'X不可逆或为了防止过拟合，我们增加lambda扰动，得到

从另一个角度来看，这相当与给我们的线性回归参数增加一个惩罚因子，这是必要的，我们数据是有干扰的，不正则的话有可能数据对于训练集拟合的特别好，但是对于新数据的预测误差很大。

4、梯度下降法

我们在上边给出了最小二乘法求解线性回归的参数theta,实际python的numpy库就是使用的这种方法。当然了，这需要我们的参数的维度不大，当维度大的时候，使用解析解就不适用了，矩阵求逆的复杂度是O（n^3），非常慢，这里讨论梯度下降算法。

上面已经说明该误差函数是凸函数，具有全局最小值

梯度下降法公式为：

对损失函数求导可以得到

这里应该是省略了前面那个sigma。使用该公式带入梯度下降原始公式可以构成梯度下降法求解线性回归的核心迭代公式。

4.1随机梯度下降法（stachastic gradient descend，SGD）

随机梯度下降法就是每次取一个样本来计算，每一次的计算量都很小，而且好处是，当数据量很多的时候，可能只用了一部分数据就能让损失达到最小。

4.2批量梯度下降法（batch gradient descend，BGD）

随机梯度下降法是逐样本一个个进行迭代的，由于数据中难免存在很多噪声，并不能保证梯度每次都是下降的，解空间的搜索是盲目进行的（工程上也是可行的，因为总体上是下降的），因此这里可以一次性对全部样本进行计算：

工程上，这里最好对误差部分求取一下均值，可以理解为梯度走动的平均值，具有稳定损失的作用。不过该方法每次迭代都是用全部的训练数据，且一次迭代一般都达不到最优解也无法判断是否最优解，计算复杂度很高。

4.3小批量梯度下降法

每次随机取若干数据来计算，避免噪声干扰问题，也避免了一次性计算量过大的问题，能够加速损失函数的下降速度。

综合来看，如果数据量不大可以用批量梯度下降法，如果数据量大则可以使用小批量梯度下降法。

5、代码实现

5.1 一元线性回归

# x的输入是n行1列，w是一行一列的权值，bias是值，noise_sigma是噪声标准差
# 目的是根据一元最小二乘法求得直线w、b，来拟合x、w_ideal、bias、noise所确定的直线
def linear_regression_one_variable(x_data, w_ideal, bias, noise_sigma):
    # 这里主要是为了保证本函数的接口与其他函数一样
    x = x_data.reshape(x_data.shape[0])
    # 计算真实的y值，加一点干扰
    ydata = x * w_ideal[0] + bias + np.random.normal(0, noise_sigma, x.shape)
    # 利用一元最小二乘法的公式计算w和b
    m = x.shape[0]
    x_mean = x.mean()
    x2_sum = np.sum(x ** 2)
    sigma_x_mean = (x.sum())**2 / m
    w = np.sum(ydata * (x - x_mean)) / (x2_sum - sigma_x_mean)
    b = np.sum(ydata - w * x) / m
    print("w,b=", w, b)
    plt.close('all')
    plt.figure('aa')
    plt.plot(xdata, ydata, 'bo')
    plt.plot(xdata, x*w+b, 'r')
    plt.show()

随机产生若干x轴的数值，调用函数输入x数值和理想状态的w和b，最后对样本数据加入一点扰动，该值是噪声的标准差

最后结果显示如下：

5.2 梯度下降法

# 利用梯度下降法做线性回归
def linear_regression_gradient_descend(x_data, w_ideal, bias, noise_sigma, learning_rate):
    ydata = x_data.dot(w_ideal)  # np.square(x_data) + noise
    noise = np.random.normal(0, noise_sigma, ydata.shape)
    ydata = ydata + bias + noise
    ydata = ydata.reshape(ydata.shape[0], 1)  # 这里需要将ydata的形状调为跟后面y的一样，n行1列的二维矩阵

    x = np.c_[x_data, np.ones(x_data.shape[0])]  # 在最后一列全部添加1，用来增加w0或者说偏置项b

    print("随机梯度下降")
    w = np.zeros(x.shape[1], dtype=np.float64)
    for _iter in range(5000):
        index = np.random.randint(0, x.shape[0], 1)[0]
        xd = x[index]
        y0 = xd.dot(w.T)
        err = (ydata[index] - y0) * xd
        w += learning_rate * err.T
        if _iter % 4999 == 0:
            loss = np.mean(np.square(ydata - x.dot(w.reshape(w.shape[0], 1))))
            print("iter ", _iter, ":", w, loss)
    print(w)

    print("批量梯度下降")
    w = np.zeros([x.shape[1], 1], dtype=np.float64)
    for _iter in range(5000):
        y0 = x.dot(w)
        err = np.mean((ydata - y0) * x, axis=0)  # 这里求了均值，公式是sum
        e2 = err.reshape(x.shape[1], 1)
        w += learning_rate * e2
        if _iter % 4999 == 0:
            loss = np.mean(np.square(ydata - x.dot(w)))
            print("iter ", _iter, ":", w.reshape(w.shape[0]), loss)

    print(w.reshape(w.shape[0]))

    plt.close('all')
    plt.figure('aa')
    plt.plot(xdata, ydata, 'bo')
    plt.plot(xdata, x.dot(w), 'r')
    plt.show()

这里增加了学习率。

（1）数据处理

输入数据xdata是n行m列的矩阵，一共有n个样本，每个样本m个特征，w是理想系数矩阵，是一个一维向量，bias是理想偏置值，noise_sigma是噪声标准差，最后是学习率。函数一开始先根据给定的理想系数和xdata计算原始样本。最后为了计算方便，在xdata的右边增加一列，值全部为1，这样将偏置并入到w中，后续回归计算就不需要单独考虑b了。比如：

（2）随机梯度下降法

每次使用一个样本参与计算，单独列出代码：

print("随机梯度下降")
w = np.zeros(x.shape[1], dtype=np.float64)
for _iter in range(5000):
    index = np.random.randint(0, x.shape[0], 1)[0]
    xd = x[index]
    y0 = xd.dot(w.T)
    err = (ydata[index] - y0) * xd
    w += learning_rate * err.T
    if _iter % 4999 == 0:
        loss = np.mean(np.square(ydata - x.dot(w.reshape(w.shape[0], 1))))
        print("iter ", _iter, ":", w, loss)
print(w)

首先定义一个向量w，长度与x的列数一样。迭代次数设置为5000，也可以设置为其他数值，迭代收敛所需次数一般随着w的维度增加而增加，也与样本数量有关。每次迭代，首先随机选取一个样本，进行y0和误差的计算，然后更新权值w。技术上没有什么难点。略去不细说。

iter 0 : [ 0.00640153 0.01087224] 4.54779075993
iter 4999 : [ 0.49841872 0.7947959 ] 0.00234022986999

[ 0.49841872 0.7947959 ]

（3）批量梯度下降法

print("批量梯度下降")
w = np.zeros([x.shape[1], 1], dtype=np.float64)
for _iter in range(5000):
    y0 = x.dot(w)
    err = np.mean((ydata - y0) * x, axis=0)  # 这里求了均值，公式是sum
    e2 = err.reshape(x.shape[1], 1)
    w += learning_rate * e2
    if _iter % 4999 == 0:
        loss = np.mean(np.square(ydata - x.dot(w)))
        print("iter ", _iter, ":", w.reshape(w.shape[0]), loss)

print(w.reshape(w.shape[0]))

由于x是一个二维向量，在矩阵运算中逐个样本计算。因此定义w为列向量，与随机梯度下降法略有不同，主要是因为参与计算的样本的维度不同，随机梯度下降法每次是一个一维的样本参与计算，而这里是一个二维向量。

求取误差后取了误差的均值，然后进行权值的更新。这里值得一提的是，err计算出来后是一个一维向量，因此还需要把它转换为与w形状一样的e2。

实际测试中，BGD的收敛速度并不比SGD快，不过BGD得到的权值更加稳定，基本上每次都收敛到相近的值，而SGD每次运行收敛到的值都可能不太一样。

小批量随机梯度下降法与批量梯度下降法差不多，主要是为了解决样本数量过多的时候，一次性计算量过大的问题，能够在平滑噪声干扰的同时加速收敛速度。

下面展示下算法在5维变量情形下回归的效果(迭代次数不变，迭代次数多一点效果会好一点)：

xdata = np.random.rand(500, 5) * 5
linear_regression_gradient_descend(xdata, [0.5, 1.0, 2.2, 0.66, 4.8], 0.8, 0.05, 0.01)

随机梯度下降
iter 0 : [ 1.02734125 0.82701757 0.64041947 0.12654751 0.55326057 0.23822491] 287.25124087
iter 4999 : [ 0.49397389 0.99176876 2.19759631 0.66308721 4.79629027 0.81560628] 0.00342948867226
[ 0.49397389 0.99176876 2.19759631 0.66308721 4.79629027 0.81560628]
批量梯度下降
iter 0 : [ 0.60107449 0.62583283 0.62943039 0.61500329 0.70028998 0.23895635] 283.22224393
iter 4999 : [ 0.49858484 0.99780882 2.19853717 0.66096431 4.79773009 0.81729313] 0.00243657014551
[ 0.49858484 0.99780882 2.19853717 0.66096431 4.79773009 0.81729313]

5.3 TensorFlow实现线性回归

tensorfow实现线性回归就很简单了，我们不需要去理会具体的优化过程，只需要写出损失函数，把它放到指定的优化器就好了.

# 输入样本集， 理想状态的参数矩阵， 偏置， 噪声标准差
# 函数会求解样本的最小二乘曲线
# x_data ：一行一个样本，是一个二维矩阵
# w_deal:一维矩阵
def linear_regression_tf(x_data, w_ideal, bias, noise_sigma):  # x_data, w_ideal, bias, noise_sigma
    ydata = x_data.dot(w_ideal)  # np.square(x_data) + noise
    noise = np.random.normal(0, noise_sigma, ydata.shape)
    ydata = ydata + bias + noise
    ydata = ydata.reshape(ydata.shape[0], 1)  # 这里需要将ydata的形状调为跟后面y的一样，n行1列的二维矩阵

    x = np.c_[x_data, np.ones(x_data.shape[0])]  # 在最后一列全部添加1，用来增加w0或者说偏置项b
    w = tf.Variable(tf.zeros([x.shape[1], 1], dtype=tf.float64))

    y = tf.matmul(x, w)
    loss = tf.reduce_mean(tf.square(ydata - y))

    optimizer = tf.train.GradientDescentOptimizer(0.1)
    train = optimizer.minimize(loss)
    init = tf.global_variables_initializer()
    with tf.Session() as sess:
        sess.run(init)
        for step in range(2000):
            sess.run(train)
            if step % 20 == 0:
                __w = np.array(sess.run(w))
                __w = __w.reshape(__w.shape[0])
                print(step, __w, sess.run(loss))
        _w = np.array(sess.run(w))
        _w = _w.reshape(_w.shape[0])
        plt.close('all')
        plt.figure('aa')
        plt.plot(x_data, ydata, 'bo')
        plt.plot(x_data, x.dot(_w), 'r')
        # plt.plot(x1, y1, 'bo')
        # plt.plot(x2, y2, 'ro')
        plt.show()

测试输入如下：

xdata = np.random.rand(500, 5)
linear_regression_tf(xdata, [0.5, 1.0, 2.2, 0.66, 4.8], 0.8, 0.05)

输出(前三个和最后三个)：