回归其实就是通过训练数据样本训练一条具有一定泛化能力，能够拟合大部分样本数据的曲线的过程。在二维坐标系中，只有x轴和y轴，将x轴视为数据集的特征，对应的y值即为训练集的结果集。通过训练得到一条直线y = ax + b的过程就称为回归。
线性回归
在实际情况下，样本数据不一定完全的分布在直线上。
通过训练得到的回归方程，如果希望预测的准确度尽可能高，则需要在训练过程中让训练样本和回归曲线的偏差尽可能小。这里定义一个衡量这种偏差的方式，称为代价，对应的代价函数(J)：
带有^的y指的是通过回归函数计算得到的预测试，而与之相减的则是实际的样本结果y。这只是仅有一个样本的情况，实际的数据集可能包含多个样本，这就是y上角标i的作用，因此代价函数改为：
m就是样本个数，有的教程中会将所有样本预测值和实际值平方之和除以2m这个只是为了后续梯度下降时求导方便，本质是没有差别的，因为无论m还是2m都是常数，只是让代价在一个常数范围内成倍数的变化。
如何使得代价函数最小？通过观察公式。特征集x中，每个xi对应一个系数，称为θi，所有θi组成的系数向量称为θ，因此预测值y等于θx=θ0x0 + θ1x1+...+θnxn，现在考虑只有一个特征的情况。代价函数即为（θ0 + θ1x1 - y）^2。
说明一下θ0，正常的回归曲线y = ax + b，可以将b理解为θ0，对应的x0为1，因此在实际运算的过程中，往往会先将测试数据特征集x最左侧加一列，加的一列值全部为1。这样在矩阵相乘的时候就可以很好的兼容θ0项。
不难发现，只有一个特征的情况下，代价函数是一个开口向上的抛物线。高中数学的知识，该抛物线的最低点对应着抛物线方程的极小值点，也就可以理解为代价函数的最小值点。我们的工作就是求出这个代价函数的极小值点所对应的θ的值，这个过程即为梯度下降。
梯度下降
假设当前的θ对应的代价函数图像上的点如图所示。
当前位置不在极小值点，即代价偏高，需要优化θ值。为了到达代价函数的极小值点，需要当前位置向下移动一定步长，如图
通过不断下降到达极小值点，每次下降的幅度是多少。这里定义每次下降的幅度为η * ▽J。η是指定的学习率，▽J就是下降的程度，可以理解为速度，速度*时间（η）也就是下降的幅度。根据导数的知识，抛物线上一点的导数就是抛物线在该点的斜率，即可以理解为物理中的速率，该点对应的导数的值也就是我们要求的▽J。梯度下降就是当前θ不断下降η * ▽J步长的过程，因此，每次更新θ的公式：θ = θ - η * ▽J。
如果下降的幅度合适，通过不断下降，我们会逐渐接近极小值点，但往往不能真正的到达极小值点，在这种情况下，我们可以将更新θ值前后的代价函数值做差，当代价低于一个很小的值时，我们也可以认为已经基本上到达极小值点。因为越接近极小值点，斜率越小，下降的幅度就越小，代价函数的变化也就越小。
基于上面的推导，模拟梯度下降过程。
import matplotlib.pyplot as plt
import numpy as np
'''模拟梯度下降曲线'''
plt_x = np.linspace(-1, 6, 141)
#模拟一条曲线 y=(x - 2.5)^2 -1
plt_y = (plt_x - 2.5) ** 2 - 1
plt.plot(plt_x, plt_y)
plt.show()
这里直接定义的一个代价函数，在-1到6的区间上均匀的取140个点，代价函数定义为y=(x - 2.5)^2 -1，并画出图像。
根据这个定义代价求函数值的函数和代价函数导数值的函数
'''代价函数的导数'''
def dJ(theta):
return 2 * (theta - 2.5)
'''代价函数'''
def J(theta):
return (theta - 2.5)**2 - 1
模拟梯度下降过程
'''梯度下降过程'''
eta = 0.1 #学习速率
epsilon = 1e-8 #当两次下降后的差值小于epsilon，则可以认为到达了极小值点
theta = 0.0
while True:
gradient = dJ(theta)
last_theta = theta
theta -= eta * gradient
if J(last_theta) - J(theta) < epsilon:
break
print(theta)
print(J(theta))
指定学习率为0.1，epsilon就是之前说的可以认为到达极小值点的阀值，θ对应theta变量，初始化为0。打印输出
2.499891109642585
-0.99999998814289
根据高中数学的知识，y=(x - 2.5)^2 -1的极小值点应该在x=2.5的位置，这里看到求得极小值点对应的θ为2.499891109642585已经很接近2.5，可以理解为下降成功，对应的代价-0.99999998814289，如果完全到达极小值点，对应的代价应该是0。
对这个下降过程图像化展示，只需要稍微改动下代码
'''绘制梯度下降曲线'''
eta = 0.1 #学习速率
epsilon = 1e-8 #当两次下降后的差值小于epsilon，则可以认为到达了极小值点
theta = 0.0
theta_history = []
while True:
gradient = dJ(theta)
last_theta = theta
theta -= eta * gradient
theta_history.append(last_theta)
if J(last_theta) - J(theta) < epsilon:
break
plt.plot(plt_x, J(plt_x))
plt.plot(np.array(theta_history), J(np.array(theta_history)), color='r', marker='+')
plt.show()
使用theta_history记录每次的θ值，获得的图像
可以看出逐渐接近极小值点。
为了方便使用，将上述过程封装函数。
'''上述过程封装函数'''
def gradient_descent(init_theta, eta, epsilon=1e-8):
theta = init_theta
theta_history = []
while True:
last_theta = theta
theta -= eta * dJ(theta)
theta_history.append(theta)
if abs((J(last_theta) - J(theta))) < epsilon:
break
return theta_history
def plot_gradient(theta_history):
plt.plot(plt_x, J(plt_x))
plt.plot(np.array(theta_history), J(np.array(theta_history)), color='r', marker='+')
plt.show()
现在讨论下学习率η，这个值的大小直接决定了学习的效果，现在将η设置为一个很小的值。
'''设置学习速率为0.01'''
eta = 0.01
theta_history = gradient_descent(0., eta)
plot_gradient(theta_history)
设置为0.01，看看下降效果
可以看出由于η变小，每步的下降幅度也变小了，'+'边的非常密集，但也最终到达极小值点。
如果η调大会怎么样，这里将η增大为0.8。
'''当学习速率偏大的情况'''
'''这个例子说明学习速率偏大会造成影响，但只要不过大，还是可以到达极小值点'''
eta = 0.8
theta_history = gradient_descent(0., eta)
plot_gradient(theta_history)
获得图像
前几次的下降非常明显的跳到代价函数图像对称轴的另一侧，但还好是一个下降的趋势，最终也到达了极小值点。这个例子说明学习速率偏大会造成影响，但只要不过大，还是可以到达极小值点。
当学习率过大的情况，将其设为1.1。
'''当学习速率过大的情况'''
eta = 1.1
theta_history = gradient_descent(0., eta)
plot_gradient(theta_history)
执行到这，系统抛出了异常OverflowError: (34, 'Result too large')，内存溢出了，对梯度下降函数做个修改，重载一下。
'''限定次数的梯度下降，避免因为学习速率过大造成无法跳出循环'''
def gradient_descent(theta, eta, n_iter, epsilon=1e-8):
theta_history = []
iter_num = 0
while iter_num < n_iter:
last_theta = theta
theta -= eta * dJ(theta)
theta_history.append(last_theta)
if abs((J(last_theta) - J(theta))) < epsilon:
break
iter_num += 1
return theta_history
n_iter参数设置一个值，防止函数陷入无线循环，当循环次数大于n_iter时即结束。
eta = 1.1
theta_history = gradient_descent(0., eta, 10, 1e-8)
plot_gradient(theta_history)
得到图像
经过n_iter次循环，得到的下降图像距离极小值点越来越远，这就是学习率过大造成的。
总结一下：学习率过小，会导致下降幅度小，训练时间增加。学习率过大，则会导致下降过程越来越偏离代价函数极小值点。
多元线性回归中的梯度下降法
上面模拟了只有一个特征的情况，真实的数据集往往有多个特征，涉及到多元线性回归。多元回归中，▽J是代价函数对每一个θi的导数组成的向量。
代价函数图像也可能不再是一元回归中的抛物线，以二元回归为例，假设代价函数z = x^2 + 2 * y^2，则其代价函数和梯度下降如图所示。
多元回归的目标还是使得代价函数最小，每个预测值yi都是由θXi相乘得到的
θ0可以理解为回归方程的常数项，X是特征，下标从1开始，但为了方便，通常都会在特征中加一个值为1的X0项，θ0 * X0 = θ0，值不变的情况下方便矩阵相乘。
推出相应的代价函数
相应的，可以得到代价函数对每个θ值的导数▽J。
当没有1/m时的推导过程。
加上1/m时的结果。
有时也会取1/2m，但和1/m差别不大。
线性回归的向量化
之前的代码通过for循环的方式，对各个属性求代价，线性回归可以将属性X和系数θ向量化，以矩阵乘法的方式进行求解。
左侧给出代价函数的导数值，先不与每个样本对应的属性Xi,j做乘积将其转化为右侧顶部的向量。
接下来分析X，对于每个样本，都含有n个特征，m个样本就可以组成m * n大小的矩阵X。对于系数θ可以是一个n * 1的矩阵，根据线性代数，m * n的矩阵和n * 1的矩阵相乘，会得到一个m * 1的矩阵，也就是回归方程得到的结果。
再和训练数据集的y做减法，两个m * 1的矩阵相减还是一个m*1的矩阵。
之后再完成与Xi,j做乘积的操作。这一步本质上还是与X矩阵相乘m * 1和m * n的矩阵没法直接相乘，但只需要转置一下，m * 1编程1 * m就可以了，当然这和图片当中把X转置变成n * m的矩阵再乘m * 1的矩阵是一个道理，都能得到n * 1的矩阵，再乘2/m，这个就是每个属性对应的▽J。代码实现：
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2/ len(X_b)
这种方式也成为批量梯度下降。
随机梯度下降
批量梯度下降每次循环中对数据集的每个样本都进行计算得到最终的▽J，这样做很准确但时间消耗比较大。随机梯度下降采用每次对样本数据集中随机的一个样本进行计算的方式，这样做运算成本降低，但下降曲线不会像批量梯度下降一样是一条平滑的通向极小值点的曲线。可能会像下面这张图一样。
因为是随机选取数据集中的点，因此不能保证每次都会向代价函数减少最快的方向（放映到图上就是沿圆的半径），甚至有可能是向代价函数增大的方向，但一般情况下，这种方式也可以找到极小值点。
每次θ递减的幅度也是η * ▽J，但这里的η和批量梯度下降不太一样。η的常见的计算公式给出。
i_iters是当前的循环次数，b是一个因子，作用是缓冲最开始循环时η递减过快的情况。如没有b的话，第一次循环η=1/1=1，第二次循环η=1/2=0.5降幅50%，对于训练模型而言太快了，加上这么一个因子会减弱η的递减幅度。
模拟随机梯度下降
import numpy as np
import matplotlib.pyplot as plt
import sklearn.datasets as dataset
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def fit_sgd(X_train, y_train, n_iter=5, t0=5, t1=50):
'''
随机梯度下降法
:param X_train:样本特征集
:param y_train:样本结果集
:param n_iter:这里的n_iter指的是遍历次数应该是整个样本个数n_iter倍
:param t0:学习率参数t0
:param t1:学习率参数t1
:return:
'''
'''代价函数的导数'''
def dJ_sgd(theta, X_i, y_i):
return X_i * (X_i.dot(theta) - y_i) * 2
'''代价函数(其实用不上)'''
def J(theta, X_i, y_i):
return (X_i.dot(theta) - y_i)**2
'''随机梯度下降'''
def sgd(theta_init, X_train, y_train, n_iter, t0, t1):
'''计算学习率'''
def learn_rate(cur_iter):
return t0 / (t1 + cur_iter)
theta = theta_init
m = len(X_train)
for i in range(n_iter):
'''打乱原本X_train、y_train的顺序'''
rand_index = np.random.permutation(m)
X_train_rand = X_train[rand_index]
y_train_rand = y_train[rand_index]
for j in range(m):
'''随机梯度下降'''
gradient = dJ_sgd(theta, X_train_rand[j], y_train_rand[j])
theta -= learn_rate(i * j) * gradient
return theta
X_train = np.hstack([np.ones((len(X_train), 1)), X_train])
theta_init = np.random.randn(X_train.shape[1])
theta = sgd(theta_init, X_train, y_train, n_iter, t0, t1)
return theta
X_train = np.hstack([np.ones((len(X_train), 1)), X_train])
theta_init = np.random.randn(X_train.shape[1])
theta = sgd(theta_init, X_train, y_train, n_iter, t0, t1)
return theta
随机梯度下降没有用到求代价函数值，达到循环次数后即完成下降。需要注意的是，批量梯度下降中，n_iters是循环的次数限制，而在随机梯度下降中，n_iters是样本集被循环的次数，这样做也是为了提升模型的准确度。如果有m大小的样本集，n_iters值为n，则一共需要循环m * n次，每次随机抽取一个样本进行梯度下降。
模拟一组数据
'''自建数据集'''
m = 100000
X_train = np.random.normal(size=m)
X_train = X_train.reshape(-1, 1)
y_train = 4.*X_train + 3. #+ np.random.normal(0, 3, size=m)
noise = np.random.normal(0, 3, size=m)
noise = noise.reshape(-1, 1)
y_train = y_train + noise
theta = fit_sgd(X_train, y_train, n_iter=2)
print(theta)
建一个100000个样本的数据集，曲线在y=4x + 3的基础上对每个元素增加一定的噪音干扰。通过之前定义的函数得到的θ输出
[3.01165213 4.03164546]
θ0是常数项，基本上是3，θ1是x的系数近似4，因此，下降结果相对准确。
在sklearn中也封装了随机梯度下降的工具类，在linear_model包下的SGDRegressor中。
