ML-梯度下降法的详细推导与代码实现
计算
对于线性回归，梯度下降法的目标就是找到一个足够好的向量\(\theta\)，使代价函数\(J(\theta) = \sum_{i=1}^{m}(\hat{y}-y_{i})^{2}\)取得最小值。线性回归的代价函数是关于\(\theta\)的多元函数。如下：
\[J(\theta) = \sum_{i=1}^{m}(\hat{y}-y^{i})^{2} = \sum_{i=1}^{m}(\theta x^{(i)}-y^{i})^{2}
\]
\[J(\theta) = \sum_{i=1}^{m}(y^{i} - \theta_{0} - \theta_{1}x^{(i)}_1 - \theta_{2}x^{(i)}_2 - ... - \theta_{n}x^{(i)}_n )^{2}
\]
对代价函数\(J(\theta)\)的每一个\(\theta\)分量求偏导数，
\[\frac{\partial J(\theta)}{\partial \theta_i}= \begin{bmatrix}
\frac{\partial J(\theta)}{\partial \theta_1} \\
\frac{\partial J(\theta)}{\partial \theta_2} \\
\vdots \\
\frac{\partial J(\theta)}{\partial \theta_n}
\end{bmatrix} = 2 \begin{bmatrix}
\sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b\theta) (-1) \\
\sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b\theta) (-X^{(i)}_1) \\
\vdots \\
\sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b\theta) (-X^{(i)}_n)
\end{bmatrix} = 2 \begin{bmatrix}
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) \\
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) (X^{(i)}_1) \\
\vdots \\
\sum_{i=1}^{m}(X^{(i)}_b\theta - y^{(i)}) (X^{(i)}_n)
\end{bmatrix}
\]
这里将J()除以m可以缩小梯度的值：
\[J(\theta) = \frac{1}{m}\sum_{i=1}^{m}(\hat{y}-y^{i})^{2}
\]
此时梯度为：
\[\frac{\partial J(\theta)}{\partial \theta_i}= \frac{2}{m} \begin{bmatrix}
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) \\
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) (X^{(i)}_1) \\
\vdots \\
\sum_{i=1}^{m}(X^{(i)}_b\theta - y^{(i)}) (X^{(i)}_n)
\end{bmatrix}
\]
- 此处也可以除以2m，因为可以方便求导数，不过，对结果没有什么影响，可以随意
- 这里的\(J(\theta)\)分为m个部分，因为有m个点。计算梯度的时候需要把m个点的坐标代进去关于\(\theta\)的表达式。这是梯度的标准求法。这个公式，是完整的梯度公式。即，在每个\(\theta\)方向求导数。得到的梯度向量就是下降最快的方向。
- 随机梯度下降法就是简化了这个\(J(\theta)\)，不需要再把m个点都代入进去，而是随机地抽取一个点，所以计算速度大为提高。
- 小批量梯度下降法就是中和了这两个算法，每次取batch个点代进去，计算梯度。BGD的梯度最后除以了m，这里的MBGD也可以除以batch。其实除和不除，都能计算出来，因为除的这个batch，在迭代的时候，会和\(\alpha\)相乘，只不过会对\(\alpha\)的取值有所影响。
以下是求梯度的函数，很简单，但是我总是记不住，这里作详细的介绍：
def dJ(theta, X_b, y):
res = np.empty(len(theta))
res[0] = np.sum(X_b.dot(theta) - y)
for i in range(1, len(theta)):
res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
return res * 2 / len(X_b)
首先上面矩阵中的每个元素其实也需要通过矩阵运算得出
\[\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) (X^{(i)}_1) = (X^{(1)}\theta - y^{(1)})X^{(1)}_1 + (X^{(2)}\theta - y^{(2)})X^{(2)}_1 + ... + (X^{(m)}\theta - y^{(m)})X^{(m)}_1
\]
上面的式子可以用用两个向量的内积表示为，此处是内积，不是矩阵乘法:
\[\begin{bmatrix}
\theta_0 + X^{(1)}_1\theta_1 + X^{(1)}_2\theta_2 + X^{(1)}_n\theta_n ... - y^{(1)} \\
\theta_0 + X^{(2)}_1\theta_1 + X^{(2)}_2\theta_2 + X^{(2)}_n\theta_n ... - y^{(2)} \\
\vdots \\
\theta_0 + X^{(m)}_1\theta_1 + X^{(m)}_2\theta_2 + X^{(m)}_n\theta_n ... - y^{(m)}
\end{bmatrix}
\cdot
\begin{bmatrix}
X^{(1)}_1 \\
X^{(2)}_1 \\
\vdots \\
X^{(m)}_1
\end{bmatrix}
\]
上面的式子又等价于：
\[(\begin{bmatrix}
1 & X^{(1)}_{1} & X^{(1)}_{2} & \cdots & X^{(1)}_{n} \\
1 & X^{(2)}_{1} & X^{(2)}_{2} & \cdots & X^{(2)}_{n} \\
\vdots & \vdots & \ddots & \vdots \\
1 & X^{(m)}_{1} & X^{(m)}_{2} & \cdots & X^{(m)}_{n}
\end{bmatrix} \cdot
\begin{bmatrix}
\theta_0 \\
\theta_1 \\
\vdots \\
\theta_n
\end{bmatrix} -
\begin{bmatrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(m)}
\end{bmatrix})\cdot
\begin{bmatrix}
X^{(1)}_1 \\
X^{(2)}_1 \\
\vdots \\
X^{(m)}_1
\end{bmatrix}
\]
最终可以表示为，\(X[:,i]\)是python中的切片的写法：
\[(X_b \theta - y)* X[:,i]
\]
因此，除了第一行以外，每一行的偏导数都可以表示为 \((X_b \theta - y)* X[:,i]\)，上面的代码中的循环，就是对除了第一行以外的所有行，每一行都进行一次计算。这里dot()是点乘，向量点乘之后就自动求和了。所以不需要再次求和。
for i in range(1, len(theta)):
res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
所以其实是两种情况:
\((X_b \theta - y)\) 第一行
\((X_b \theta - y)* X[:,i]\) 之后的所有行
实现
思路：先准备数据集X和y ==> 使用交叉验证选出训练样本 ==> 回归
回归的过程: 将转换为X_b(就是加上X0) ==> 设置初始化 \(\theta\) ==> 进行梯度下降(计算各个方向的偏导数-> 循环计算theta = theta - eta * dj)
class LinearRegression2:
def __init__(self):
self.coef_ = None # 表示参数，theta_[1:]
self.intercept_ = None # 表示截距 ==>theta[0]
self._thera = None # 表示完整的theta==> theta[:]
def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
"""使用梯度下降法寻找最小的代价函数"""
# 格式化X和theta，加上x0 和 theta0
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
# 调用循环的梯度下降
self._thera = self.gradient_descent(X_b, y_train, initial_theta, eta=eta, n_iters=n_iters)
self.intercept_ = self._thera[0]
self.coef_ = self._thera[1:]
return self
def gradient_descent(self, X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
theta = initial_theta
i = 0
while i < n_iters:
i += 1
lastTheta = theta # 记录上一个参数向量
dj = self.dj(theta, X_b, y)
theta = theta - eta * dj
if np.absolute(self.J(lastTheta, X_b, y) - self.J(theta, X_b, y)) < epsilon:
break
return theta
def dj(self, theta, X_b, y):
"""计算代价函数的偏导数"""
# res 是 长度为len的一维数组
res = np.zeros((len(theta))) ##### zeros(5) == zeros((5,)) != zeros((5,1))
res[0] = np.sum(X_b.dot(theta) - y) # 需要将向量求和。没有办法，这里只能手动分开求解。因为这里X是从列方向分割，没有X0
for i in range(1, len(theta)):
res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
return res * 2 / len(X_b)
def J(self, theta, X_b, y):
"""计算代价函数"""
return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
调用效果如下：
np.random.seed(666)
x = 2 * np.random.random(size=100)
y = x * 3. + 4. + np.random.normal(size=100)
X = x.reshape(-1, 1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_ratio=0.2, seed=666)
line = LinearRegression2()
line.fit_gd(X_train, y_train)
###
4.085675667692203
[ 2.97732994]
使用向量实现
博客没法显示全公式。。就用图片凑合吧.\(X_b\theta -y\)是列向量
即，最后的公式为：
\[J(\theta) = \frac{2}{m}X^T_b \cdot (X_b\theta-y)
\]
使用python的numpy实现：
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)
其中，\(X_b\)是在原始数据的每一列前面加了一列后的矩阵，y是训练数据中给定的\(y\_train\)。
小结
- 不管是用向量计算还是用手动计算。都需要先把训练数据\(X\)转换为\(X_b\)
- 一般来说，\(X\)都是ndarray类型的二维数组，使用
hstack([ist1, list2])函数可以将两个数组摊在一起。将\(X\)转换为\(X_b\)的代码为：X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
- 大体思路就是:
- 调用
fit()方法拟合，该方法中产生\(X_b\)，并初始化\(\theta\)。最后使调用梯度下降方法gradient_descent()来找到最优的\(\theta\) gradient_descent()循环调用dj()计算对\(\theta\)的偏导数,并每一次都对\(\theta\)的值进行更新。gradient_descent()中会调用J()来判断梯度的增量是否已经足够小。
- 调用
数据标准化
在这里使用波士顿房价数据测试，发现一些很有意思的东西，如下。为什么会这样？，因为数据没有归一化，在数据集中存在一下特别大和特别小的数字，导致在计算梯度时，可能会导致梯度的跨度太大，而无法收敛。也有可能时计算式出现了过大的数inf。
linear.fit_gd(X, y, n_iters=1e4, eta=0.001)
print(linear.coef_)
###
[nan nan nan nan nan nan nan nan nan nan nan nan nan]
只要将学习系数设定的足够小，就不会报错。
linear.fit_gd(X, y, n_iters=1e4, eta=0.0000001)
使用sklearn的StandardScaler归一化数据。注意归一化需要将所有的X都归一化，包括用来训练的x_train和用来测试的x_test
from sklearn.preprocessing import StandardScaler
standardScaler = StandardScaler()
standardScaler.fit(X) # 导入数据
X_standard = standardScaler.transform(X) # 转换数据
所以，到目前为止，要预测一个波士顿房价数据的完整过程是这样的：
from sklearn import datasets
# 加载数据集
X = datasets.load_boston().data
y = datasets.load_boston().target
# 剔除噪音
X = X[y < 50]
y = y[y < 50]
# 数据归一化处理
from sklearn.preprocessing import StandardScaler
standardScaler = StandardScaler()
standardScaler.fit(X)
X_standard = standardScaler.transform(X)
# 交叉验证
X_train_standard, X_test_standard, y_train, y_test = train_test_split(X_standard, y, test_ratio=0.2, seed=666)
# 回归，拟合
linear = LinearRegression2()
linear.fit_gd(X_train_standard, y_train, n_iters=3e5)
# 使用score计算拟合的效果
print(linear.score(X_test_standard, y_test))
随机梯度下降法
原本的梯度下降法：
\[\frac{\partial J(\theta)}{\partial \theta_i}= \frac{2}{m} \begin{bmatrix}
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) \\
\sum_{i=1}^{m}(X^{(i)}_b\theta -y^{(i)}) (X^{(i)}_1) \\
\vdots \\
\sum_{i=1}^{m}(X^{(i)}_b\theta - y^{(i)}) (X^{(i)}_n)
\end{bmatrix}
\]
随机梯度下降法：
\[\text { 2. }\left(\begin{array}{c}{\left(X_{b}^{(i)} \theta-y^{(i)}\right) \cdot X_{0}^{(i)}} \\
{\left(X_{b}^{(i)} \theta-y^{(i)}\right) \cdot X_{1}^{(i)}} \\
{\left(X_{b}^{(i)} \theta-y^{(i)}\right) \cdot X_{2}^{(i)}} \\
{\cdots} \\
{\left(X_{b}^{(i)} \theta-y^{(i)}\right) \cdot X_{n}^{(i)}}\end{array}\right)=2 \cdot\left(X_{b}^{(i)}\right)^{T} \cdot\left(X_{b}^{(i)} \theta-y^{(i)}\right)
\]
原本的计算偏导数的函数用这个替代。里只有一个样本哦，所有的误差和X都是一个同一个样本的。而原来的函数，每计算一次偏导数都要对n个特征，每个循环m次，共mn次计算。计算量确实是大大提高，但是，不能保证每次都找到最好的梯度。使用随机产生的梯度。相比较原来的函数，去掉了步长。
学习率：
\[\eta=\frac{t_{0}}{i_{-} i t e r s+t_{1}}
\]
- 去掉了参数中的步长
- n_iters 表示遍历所有样本的轮数
- 学习率需要不断缩小，因为，算出来的梯度不能保证是越来越小的（因为随机）
class StochasticDescent:
"""随机梯度下降法"""
def __init__(self):
self.coef_ = None # 表示参数，theta_[1:]
self.intercept_ = None # 表示截距 ==>theta[0]
self._thera = None # 表示完整的theta==> theta[:]
def fit_SGD(self, X_train, y_train, eta=0.01, n_iters=1):
"""使用梯度下降法寻找最小的代价函数"""
# 格式化X和theta，加上x0 和 theta0
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
# 调用循环的梯度下降
self._thera = self.sgd(X_b, y_train, initial_theta, n_iters, 5, 50)
self.intercept_ = self._thera[0]
self.coef_ = self._thera[1:]
return self
def sgd(self, X_b, y_train, initial_theta, n_iters, t0=5, t1=50):
"""随机梯度下降, niters是轮数"""
def get_study_rate(i_iters):
"""把学习率和迭代次数联系起来"""
return t0 / (t1 + i_iters)
def dJ_sgd(X_b_i, theta, y):
"""计算随机一个元素的梯度"""
return 2 * X_b_i.T.dot(X_b_i.dot(theta) - y)
m = len(X_b)
theta = initial_theta
for n in range(int(n_iters)):
indexes = np.random.permutation(m)
X_b_new = X_b[indexes, :]
y_new = y[indexes]
for i in range(m):
sgd = dJ_sgd(X_b_new[i], theta, y_new[i])
theta = theta - get_study_rate(m * n + i) * sgd
return theta
sklearn 中的随机梯度下降
SGDRegressor位于线性模型下，SGDRegressor(max_iter=50)可以传入参数 max_iter指定迭代的次数。
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50)
sgd_reg.fit(X, y)
调试梯度下降法
对每个\(\theta\)分量，都要微积分，最后，用算出来的微积分向量，作为梯度。
- 计算量巨大，调试成功后，要关闭
- 分子的两个顺序不能错，我把顺序搞错了，结果导数算反了。。。
\[\frac{d J}{d \theta}=\frac{J(\theta+\varepsilon)-J(\theta-\varepsilon)}{2 \varepsilon}
\]
关键代码在这：
def dj_debug(self, X_b, theta, y, epsilon=0.001):
res = np.zeros((len(theta), ))
for i in range(len(theta)):
theta_1 = theta.copy()
theta_2 = theta.copy()
theta_1[i] += epsilon
theta_2[i] -= epsilon
res[i] = (self.J(theta_1, X_b, y) - self.J(theta_2, X_b, y)) / (2 * epsilon)
return res
完整代码：
class aDebugGradient:
"""随机梯度下降法"""
def __init__(self):
self.coef_ = None # 表示参数，theta_[1:]
self.intercept_ = None # 表示截距 ==>theta[0]
self._thera = None # 表示完整的theta==> theta[:]
def fit_debug(self, X_train, y_train, eta=0.01, n_iters=1e4):
"""使用梯度下降法寻找最小的代价函数"""
# 格式化X和theta，加上x0 和 theta0
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
# 调用循环的梯度下降
self._thera = self.gradient_Descent(self.dj_debug, X_b, y_train, initial_theta, eta, n_iters)
self.intercept_ = self._thera[0]
self.coef_ = self._thera[1:]
return self
def J(self, theta, X_b, y):
return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
def dj_debug(self, X_b, theta, y, epsilon=0.001):
res = np.zeros((len(theta), ))
for i in range(len(theta)):
theta_1 = theta.copy()
theta_2 = theta.copy()
theta_1[i] += epsilon
theta_2[i] -= epsilon
res[i] = (self.J(theta_1, X_b, y) - self.J(theta_2, X_b, y)) / (2 * epsilon)
return res
def dJ_math(self, X_b, theta, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)
def gradient_Descent(self, dj, X_b, y_train, initial_theta, eta, n_iters, epsilon=1e-8):
"""梯度下降"""
theta = initial_theta
for i in range(int(n_iters)):
gradient = dj(X_b, theta, y_train) # 计算得到梯度
last_theta = theta
theta = theta - eta * gradient
if np.absolute(self.J(theta, X_b, y_train) - self.J(last_theta, X_b, y_train)) < epsilon:
break
return theta
最后，，，这垃圾公式真是烦死人！！！
