1 线性目标的梯度优化
损失函数:
算法1 : 批量梯度下降BGD
每次迭代使用所有样本来对参数进行更新。
损失函数:
代数形式:
矩阵形式:
更新:
代数形式伪代码:
矩阵形式伪代码:
算法2:随机梯度下降SGD
每次迭代使用一个样本来对参数进行更新。
一个样本的损失函数:
代数形式伪代码:
矩阵形式伪代码
算法3:小批量梯度下降法
每次迭代使用n_batch个样本对参数进行更新
代数形式伪代码:
矩阵形式伪代码
1. 线性目标的sgd优化,写出伪代码。mini_batch怎么更新?
2. SGD并行实现?
=============================
# 梯度优化代码
import numpy as np import pandas as pdfrom sklearn.model_selection import train_test_split from sklearn.linear_model import LinearRegression #基准模型 data_df = pd.read_csv('kc_house_data.csv') train_data, test_data = train_test_split( data_df, train_size=0.8, random_state=1234 ) # sklearn模型 X_train = train_data['sqft_living'].values.reshape(-1,1) y_train = train_data['price'].values.reshape(-1,1) X_test = test_data['sqft_living'].values.reshape(-1,1) y_test = test_data['price'].values.reshape(-1,1) model = LinearRegression(normalize=True).fit(X_train, y_train ) y_pred = model.predict(X_test) print('intercept:',model.intercept_) print('coel:', model.coef_) print('MSE: %.2f' % np.mean((y_pred - y_test) ** 2)) # 模型 def linear(X,w): X = np.hstack((np.ones((X.shape[0],1)), X)) y = np.dot(X,w) return y # 损失函数 def cost_function( X, y, w ): m = X.shape[0] e = np.dot(X,w) - y J = 0.5/m * np.dot(e.T, e) return J[0] # 标准化 def normalize( X ): mu = np.mean(X,axis =0) sigma = np.std(X,axis = 0) X = (X-mu)/sigma return X,mu,sigma # 均方误差 def MSE( X, y, w ): return np.mean((linear(X,w) - y) ** 2) #----------------------------- # 批量梯度下降法 def batch_gradient_descent( X, y, eta=0.1, eps=1e-4, num_iters=2000 ): converged = False J = np.zeros( num_iters ) m = X.shape[0] i_iters = 0 w = np.zeros((X.shape[1],1)) while not converged and i_iters < num_iters: grad_w = np.dot( X.T, X.dot(w) - y ) w -= eta*grad_w/m dist_grad = np.sqrt(np.dot(grad_w.T,grad_w)) J[i_iters] = cost_function( X, y, w ) print('Iteration:', i_iters+1, ': J(w):', J[i_iters], "||J(w)||:", dist_grad) if dist_grad < eps: cur= True i_iters += 1 return w X, mu, sigma = normalize( X_train ) X = np.hstack( (np.ones((X_train.shape[0],1)), X) ) w = batch_gradient_descent( X, y, eta = 0.1, eps=1e-4 ) print ('Intercept:', w[0] - np.dot(mu/sigma,w[1:])) print( 'Slope :', w[1:]/sigma[:,np.newaxis]) print ("MSE : %.2f" % MSE( (X_test-mu)/sigma, y_test, w )) #---------------------------------------------- # 随机梯度下降法 def stochastic_gradient_descent( X, y, eta=0.1, eps=1e-4, num_iters=2000 ): cur = False J = np.zeros( num_iters ) n_batch = 1 ## i_iters = 0 w = np.zeros((X.shape[1],1)) while not cur and i_iters < num_iters: batch = np.random.choice(X.shape[0], n_batch) X_batch, y_batch = X[batch], y[batch] grad_w = np.dot( X_batch.T, X_batch.dot(w) - y_batch ) w -= eta*grad_w/n_batch # dist_grad = np.sqrt(np.dot(grad_w.T,grad_w)) J[i_iters] = cost_function( X_batch, y_batch, w ) print('Iteration:', i_iters+1, ': J(w):', J[i_iters], "||J(w)||:", dist_grad) if dist_grad < eps: cur = True i_iters += 1 return w X, mu, sigma = normalize( X_train ) X = np.hstack( (np.ones((X_train.shape[0],1)), X) ) w = stochastic_gradient_descent( X, y, eta = 0.1, eps=1e-4 ) print ('Intercept:', w[0] - np.dot(mu/sigma,w[1:])) print( 'Slope :', w[1:]/sigma[:,np.newaxis]) print ("MSE : %.2f" % MSE( (X_test-mu)/sigma, y_test, w )) # --------------------------------------- # 小批量随机梯度下降法 def mini_gradient_descent( X, y, eta=0.1, eps=1e-4, num_iters=2000 ): cur = False J = np.zeros( num_iters ) n_batch = 5000 ## i_iters = 0 w = np.zeros((X.shape[1],1)) while not cur and i_iters < num_iters: batch = np.random.choice(X.shape[0], n_batch) X_batch, y_batch = X[batch], y[batch] grad_w = np.dot( X_batch.T, X_batch.dot(w) - y_batch ) w -= eta*grad_w/n_batch # n_batch为所选择的批量个数 dist_grad = np.sqrt(np.dot(grad_w.T,grad_w)) J[i_iters] = cost_function( X_batch, y_batch, w ) print('Iteration:', i_iters+1, ': J(w):', J[i_iters], "||J(w)||:", dist_grad) if dist_grad < eps: cur = True i_iters += 1 return w X, mu, sigma = normalize( X_train ) X = np.hstack( (np.ones((X_train.shape[0],1)), X) ) w = mini_gradient_descent( X, y, eta = 0.1, eps=1e-4 ) print ('Intercept:', w[0] - np.dot(mu/sigma,w[1:])) print( 'Slope :', w[1:]/sigma[:,np.newaxis]) print ("MSE : %.2f" % MSE( (X_test-mu)/sigma, y_test, w ))
代码写的再多,不如自己尝试一遍,很多东西都是自己操作了才学得比较快。这里的数据以及大部分知识点多是从公众号学习,这里就不具体说明。很多东西都是要自己找到的才会真正舍得去学。数据文件在网盘,http://pan.baidu.com/share/link?shareid=1375610667&uk=209896182。可以自己尝试下载数据练习,里面也有各种jupyter,还有岭回归和拉索回归的,感兴趣可以自行去学习。
2 LR回归原理推导
假设数据服从伯努利分布,通过极大似然函数的方法,用梯度下降来求解参数,来达到将数据分类,一般是二分类。
2.1 对率回归推导
1. sigmoid函数推导
逻辑回归解决的是二分类问题,首先通过训练模型来学习x和y的映射关系,有了映射关系就容易对x和y进行预测分类。如何定义这种映射关系,这里假定y的取值为0或者1(也可以是其他类),两种可能,自然而然服从伯努利分布,这种映射关系就可以通过条件概率p(y|x)来表示。p(y=1|x)和p(y=0|x),对于二分类就可设定一个阈值,0.5。
如何通过条件概率来描述x和y之间的关系呢?对率回归是线性模型,但是无法直接表示成,那么就可以通过广义线性模型来实现。一般广义线性模型的一般形式,
,其中g(y)为单调函数,g(y)就可被作为联系函数。(为什么选择sigmoid函数,一般先考虑的是单位阶跃函数将任意的实数转换为0/1的概率值,用它来当成联系函数来判断属于哪个类。但是单位阶跃函数的一个缺点就是在零点不连续也不单调,而联系函数需要单调连续。)而sigmoid能够将(-无穷,+无穷)的值域映射到(0,1),这样就可以得到合理的概率值,而且单调连续,可以作为联系函数。
这样就得到了,将
代入,得到
用条件概率表示就是
得到
2. 广义线性模型推导
指数族分布:是一类分布的总称,它的概率密度函数一般形式是:
其中,称为该分布的自然参数;T(y)为充分统计量,通常为y本身;
为配分函数,保证概率表达式加和为1,保证式子是合格的概率密度函数;b(y)是关于随机变量y的函数。常见伯努利和正态均为指数族分布。
证明:伯努利分布是指数族分布?
化成指数族的一般形式:
对应指数族分布的一般形式,
广义线性模型三假设:
1. 给定x的条件下,假设随机变量y服从指数族分布
2. 给定x条件下,目标是得到一个模型h(x)能预测出T(y)的期望值。
3. 假设该指数族分布中的自然参数和x呈线性关系,即
满足这三条假设即为广义线性模型。
对数几率回归是在对一个二分类进行建模,并假设被建模的随机变量y取值为0或1,可以很自然地假设y服从伯努利分布。如果想要构建一个线性模型来预测在给定x的条件下y取值的话,可以考虑使用广义线性模型来建模。
假设1:伯努利分布服从指数族分布。
假设2:,得到
假设3:
2.2 损失函数推导
损失函数:
由式子,
假设数据都是独立的,由伯努利分布,得到一个样本发生的条件概率表示为:
进而得到它们一起发生的概率(似然函数):
取对数,取负,得到损失函数
为什么用似然函数?
我们的目标是预测为某一类的出现概率最大,每个样本预测都要得到最大的概率。前面我们得到了一个样本的条件概率,极大似然估计就会将所有的样本考虑进去,来使得观测到的样本的出现概率最大,所以有了累乘形式的似然函数。 (为什么累乘?)基于条件独立的假设,总的条件概率就可以表示为每一个条件概率的累乘形式。
为什么用极大似然函数(交叉熵)作为损失函数而不用最小二乘(欧氏距离,均方损失)?
均方损失:假设误差是正态分布,适合线性的输出(回归问题),特点是对于与真实值差别越大,惩罚力度越大,不适合分类问题。
交叉熵损失:假设误差是伯努利分布,可以视为预测概率分布与真实概率分布的相似程度。多应用在分类的问题。
均方误差对参数的偏导的结果都乘以了 sigmoid的导数,而sigmoid导数在其变量值很大或很小时趋近于0,所以均方误差的偏导可能接近于0。
而参数更新公式: 参数 = 参数 – 学习率 * 损失函数对参数的偏导
在偏导很小时,参数更新速度会变得很慢,而在接近于0时,参数几乎不更新,出现梯度消失的情况。反观交叉熵对参数的偏导就没有sigmoid导数,所以不存在梯度消失的问题。
均方损失:
使用梯度下降法的条件损失函数时凸函数。而对最小二乘的损失函数求导,
可以知道,J(w)对w不是凸函数,不能用代价函数。
交叉熵:
当越接近于1,越接近于0。预测值与真实值完全相同,其损失函数为0。
当越接近于0,越接近于。
当越接近于1,越接近于。
当越接近于0,越接近于0。预测值与真实值完全相同,其损失函数为0。
写成矩阵形式:
2.3 梯度下降法
代数形式伪代码:
矩阵形式伪代码:
==================================
# LR代码
import numpy as np #随机生成样本数据,二分类问题,每类样本生成5000个 np.random.seed(43) num = 5000 x1 = np.random.multivariate_normal([0,0],[[1,0.75],[0.75,1]],num) #d多元正态分布矩阵,[0,0]为均值,然后是正定对称半正定矩阵。 x2 = np.random.multivariate_normal([1, 4], [[1, .75],[.75, 1]], num) X = np.vstack((x1, x2)).astype(np.float32) X = np.hstack( (np.ones((X.shape[0],1)), X) ) y = np.hstack((np.zeros(num),np.ones(num))) # 模型 def sigmoid(x): return 1 / (1 + np.exp(-x)) # 损失函数 def log_likelihood(X, y, w): # 首先按照标签来提取正样本和负样本的下标 pos, neg = np.where(y==1), np.where(y==0) pos_sum = np.sum(np.log(sigmoid(np.dot(X[pos], w)))) neg_sum = np.sum(np.log(1-sigmoid(np.dot(X[neg], w)))) return -(pos_sum + neg_sum) # 返回cross entropy loss # 实现逻辑回归模型 def logistic_regression(X, y, num_steps, eta): m = X.shape[1] w = np.zeros(X.shape[1]) for step in range(num_steps): error = sigmoid(np.dot(X,w)) - y grad_w = np.matmul(X.T, error) w = w - eta/m * grad_w # 每隔一段时间,计算一下损失函数,看看有没有变化 # 正常情况下, 它会慢慢变小,最后收敛 if step % 10000 == 0: print (log_likelihood(X, y, w)) return w w = logistic_regression(X, y, num_steps = 100000, eta = 5e-5) print ("参数w, b分别为: ", w[1:],w[0])# sklearn from sklearn.linear_model import LogisticRegression # C设置一个很大的值,意味着不想加入正则项 clf = LogisticRegression(fit_intercept=True, C = 1e15) clf.fit(X, y) print ("(sklearn)逻辑回归的参数w, b分别为: ", clf.coef_, clf.intercept_, )
=============================================================================================
写在后,这是第一次写博客,平常喜欢在word或txt记录,也习惯了word的公式,就不打latex了,公式直接截图了。发出来的初衷有两个,一个是怕万一本地的文档像上个电脑突然去世,里面的资料都丢失了。还一个也是希望,对上述公式或者知识点是以个人的见解来总结或者直接拿过来,希望能够看到大家的看法,以及纠正,在这也对直接拿过来的这些原出处没办法做引用说明说声抱歉。
也是在自学的路上,更新会比较慢,这也是对自己的一种鞭挞吧。
后续,还会有想对牛顿法,拟牛顿法进行补充,在能够认为掌握之后。
愿正在阅读的你,以及仍热爱的你,共勉。