machinelearning之梯度下降（bgdsgd）

问题 简单来说，有一堆实数数据，数据的格式如下： x1,x2,x3,?,xn,y 所有的这些数据称为训练集，其中x称为feature，y称为target。 现在又有一些数据： x1,x2,x3,?,xn 需要做的是根据这些x的值，推测出y的值。 解决方法 Overdetermined Equations 假设y是x的

问题

简单来说，有一堆实数数据，数据的格式如下：

x1,x2,x3,?,xn,y

所有的这些数据称为训练集，其中x称为feature，y称为target。

现在又有一些数据：

x1,x2,x3,?,xn

需要做的是根据这些x的值，推测出y的值。

解决方法

Overdetermined Equations

假设y是x的线性函数（顺便说一句lr中的linear是对于θ而言的，并非针对x），表达为公式为：

y=θ0x0+θ1x1+θ2x2+?+θnxn

其中x0为截距(intercept term)，其值恒为1。

最容易想到的方法，可以把所有训练集的数据代入这个公式，得到方程组：

???????????????y(1)=θ0x(1)0+θ1x(1)1+θ2x(1)2+?+θnx(1)ny(2)=θ0x(2)0+θ1x(2)1+θ2x(2)2+?+θnx(2)n?y(m)=θ0x(m)0+θ1x(m)1+θ2x(m)2+?+θnx(m)n

这个方程组有m个方程，n+1个未知数，实际问题中通常是训练集的个数大于feature个数，也就是说m > n+1，这种情况下的方程组称为超定方程组，是不能直接求解的。当然可以像当年欧拉和拉普拉斯最初解决天文计算问题一样(here)，把m个方程组分成n+1组，然后每一组合并成一个方程，得到n+1个方程后再求解。不过问题是怎么分成n+1组，这个很是adhoc的。

Cost Function

机器学习上解决这个问题的方法是定义一个损失函数：

J(θ)=12∑i=1m(hθ(x(i))?y(i))2

然后选择适当的θ，使得J(θ)最小。

BatchGradient Descent

这个最小化的算法在机器学习中称为梯度下降：

•随机初始化一组θ值；

•朝着减少cost function的方向，不断更新θ值，直到收敛。更新公式为： θj:=θj?α?J(θ)?θj

其中α为学习速率(learning rate)。

Gradient Descent推导

假设训练集中只有一个数据，?J(θ)?θj计算如下：

?J(θ)?θj=?(12(hθ(x)?y)2)?θj=2?12(hθ(x)?y)??(hθ(x)?y)?θj=(hθ(x)?y)??(hθ(x)?y)?θj=(hθ(x)?y)??(∑ni=0θixi?y)?θj=(hθ(x)?y)xj

代入更新公式：

θj=θj?α(hθ(x)?y)xj=θj+α(y?hθ(x))xj

对于有m个数据集的情况可以得到如下公式：

θj:=θj+α∑i=1m(y(i)?hθ(x(i)))x(i)j

Gradient Descent直观解释

J(θ)是一个关于θ的多元函数，高等数学的知识说，J(θ)在点P(θ0,θ1,?,θn)延梯度方向上升最快。现在要最小化 J(θ)，为了让J(θ)尽快收敛，就在更新θ时减去其在P点的梯度。

在最终推导出的更新公式中，可以得出以下直观结论：如果遇到一个数据使得(y?hθ(x))比较小，这时候θ的更新也会很小，这也符合直观感觉。当一个数据使得差值比较大时，θ的更新也会比较大。

Stochastic Gradient Descent

以上的讨论的算法叫batch gradient descent，batch指的是，每次更新θ的时候都需要所有的数据集。这个算法有两个缺陷：

◦数据集很大时，训练过程计算量太大；

◦需要得到所有的数据才能开始训练；

比如一个场景下，我们训练了一个lr模型，应用于线上环境，当这个模型跑在线上的时候我们会收集更多的数据。但是上面两个问题使得我们不能及时更新模型，而这正是随机梯度下降要解决的问题。

在之前的推导过程中已经给出了sgd的更新公式，只是没有指出，现正式提出sgd的更新公式：

loop for every (x, y) in training set until convergence:

θj:=θj+α(y?hθ(x))xj
与bgd唯一的区别是，无论数据集有多少，每次迭代都只用一个数据。这样当有新的数据时，直接通过上式更新θ，这就是所谓的online learning。又因为每次更新都只用到一个数据，所以可以显著减少计算量。

 
批量梯度下降是一种对参数的update进行累积，然后批量更新的一种方式。用于在已知整个训练集时的一种训练方式，但对于大规模数据并不合适。

随机梯度下降是一种对参数随着样本训练，一个一个的及时update的方式。常用于大规模训练集，当往往容易收敛到局部最优解。

说明:因为最小二乘问题是一个求凸函数极值的问题，它只有一个最优解，没有所谓的局部最优，所以在这个问题上完全可以大用梯度下降来解

Mini-batch gradient
它还是采用了batch的思路，也就是所有样本一起更新。和batch不同的是mini，在求解方向的时候选择了一部分样本一起更新，这样就减少了计算量，同时它又不像SGD那样极端只使用一个样本，所以保证了方向的精确性。一句话总结就是，mini-batch是一个位于BGD和SGD之间的算法，精度比BGD低，比SGD高，速度比BGD快，比SGD慢（这个结论只是单从公式上分析，没有实证）。
看下面的迭代公式，则是10个一组进行更新。

附：

梯度gradient

http://zh.wikipedia.org/wiki/%E6%A2%AF%E5%BA%A6

在标量场f中的一点处存在一个矢量G，该矢量方向为f在该点处变化率最大的方向，其模也等于这个最大变化率的数值，则矢量G称为标量场f的梯度。

在向量微积分中，标量场的梯度是一个向量场。

标量场中某一点上的梯度指向标量场增长最快的方向，梯度的长度是这个最大的变化率。

更严格的说，从欧氏空间Rn到R的函数的梯度是在Rn某一点最佳的线性近似。在这个意义上，梯度是雅戈比矩阵的一个特殊情况。

在单变量的实值函数的情况，梯度只是导数，或者，对于一个线性函数，也就是线的斜率。

梯度一词有时用于斜度，也就是一个曲面沿着给定方向的倾斜程度。

一个标量函数的梯度记为： 或 ， 其中（nabla）表示矢量微分算子。

梯度下降法

http://zh.wikipedia.org/wiki/%E6%A2%AF%E5%BA%A6%E4%B8%8B%E9%99%8D%E6%B3%95

梯度下降法，基于这样的观察：

如果实值函数  在点  处可微且有定义，那么函数 在  点沿着梯度相反的方向  下降最快。因而，如果

对于  为一个够小数值时成立，那么 。

是向量。

考虑到这一点，我们可以从函数  的局部极小值的初始估计  出发，并考虑如下序列  使得

因此可得到

如果顺利的话序列  收敛到期望的极值。注意每次迭代步长  可以改变。

梯度下降法的缺点是：

•靠近极小值时速度减慢。

•直线搜索可能会产生一些问题。

•可能会'之字型'地下降。

随机梯度下降法,也叫增量梯度下降

由于梯度下降法收敛速度慢，而随机梯度下降法会快很多

根据某个单独样例的误差增量计算权值更新，得到近似的梯度下降搜索（随机取一个样例）

可以看作为每个单独的训练样例定义不同的误差函数

在迭代所有训练样例时，这些权值更新的序列给出了对于原来误差函数的梯度下降的一个合理近似

通过使下降速率的值足够小，可以使随机梯度下降以任意程度接近于真实梯度下降

?标准梯度下降和随机梯度下降之间的关键区别

标准梯度下降是在权值更新前对所有样例汇总误差，而随机梯度下降的权值是通过考查某个训练样例来更新的

在标准梯度下降中，权值更新的每一步对多个样例求和，需要更多的计算

标准梯度下降，由于使用真正的梯度，标准梯度下降对于每一次权值更新经常使用比随机梯度下降大的步长

如果标准误差曲面有多个局部极小值，随机梯度下降有时可能避免陷入这些局部极小值中

sgd、bgd的Python实现
#coding=gbk
'''
Created on Apr 12, 2014

@author: pipi
'''
import numpy as np

def bgd(feature,target,alpha = 0.001,iterateTimes = 200):
'... batch gradient descent ...'
theta = np.zeros(feature.shape[1])
for it in range(iterateTimes): #max iteratetimes is 200
for i in range(feature.shape[0]): #for each sample
error = target[i] - sum(feature[i]*theta)
theta += alpha*error*feature[i]

predict = [sum(theta*sample) for sample in feature]
mse = sum((predict - target)**2)/feature.shape[0]
print 'bgd_mse : ',mse
return theta

def sgd(feature,target,alpha = 0.001,iterateTimes = 101000):#101000
'... stochastic gradient descent ...'
theta = np.zeros(feature.shape[1])#num of theta = num of feature atrribute
for it in range(iterateTimes): #max iteratetimes is 200
i = it%feature.shape[0]
error = target[i] - sum(feature[i]*theta)#对应元素相乘，都是行array
theta += alpha*error*feature[i]

predict = [sum(theta*sample) for sample in feature]
mse = sum((predict - target)**2)/feature.shape[0]
if(mse break
print 'sgd_mse : ',mse

return theta

def normalizer(feature):
'normalization of feature'
mean_j = np.mean(feature,axis = 0)
for j in range(1,feature.shape[1]):
feature[:,j] = (feature[:,j] - mean_j[j])/std_j[j]
return feature

'''
Created on Apr 12, 2014

@author: pipi
'''
import re
import numpy as np

def loadData(filename):
feature = list()
target = list()
f = open(filename,'rb')
for line in f:
sample = re.split('\s+',line.strip())
feature.append([1] + sample[0:-1])#construct x0 = 1
target.append(sample[-1])
return np.array(feature,np.float),np.array(target,np.float)

代码中使用的数据集可以从http://download.csdn.net/detail/pipisorry/7192349下载

代码中normalize函数用于对feature进行归一化处理，可以尝试一下去掉normalize过程，对于这个数据集会得出很出乎意料的结果。
 
可能存在的改进

1）样本可靠度，特征完备性的验证

例如可能存在一些outlier，这种outlier可能是测量误差，也有可能是未考虑样本特征，例如有一件衣服色彩评分1分，料子1分，确可以卖到10000万元，原来是上面有一个姚明的签名，这个特征没有考虑，所以出现了训练的误差，识别样本中outlier产生的原因。

2）批量梯度下降方法的改进

并行执行批量梯度下降

3）随机梯度下降方法的改进

找到一个合适的训练路径（学习顺序），去最大可能的找到全局最优解

4）假设合理性的检验

H（X)是否合理的检验

5）维度放大

维度放大和过拟合问题，维度过大对训练集拟合会改善，对测试集的适用性会变差，如果找到合理的方法？

概率解释

在以上的讨论中，得出y与x的关系是线性假设，使用梯度下降也可以从高数中得到依据，唯有损失函数好像是拍脑袋想出来的。有那么多的函数可以用，为什么单选择了一个二次式做为损失函数。其实这里选择二次函数是有其理论基础的。

y与x满足以下公式：

y(i)=θTx(i)+ε(i)

其中ε(i)称为误差，可能由两个原因产生：

◦feature选择的不合适；

◦随机噪声；

又假设ε(i)独立同分布，且满足均值为0，方差为σ2的高斯分布，即：

p(ε(i))=12π??√σe?(ε(i))22σ2

也就是：

p(y(i)|x(i);θ)=12π??√σe?(y(i)?θTx(i))22σ2

以上是一个关于y, X的公式，可以定义一个似然函数，形式如同上式，但是似然函数是关于θ的公式：

L(θ)=L(θ;X,y)=p(y|X;θ)

根据之前ε(i)的独立性假设，L(θ)可以记做

L(θ)=∏i=1mp(y(i)|x(i);θ)=∏i=1m12π??√σe?(y(i)?θTx(i))22σ2

现在已经观察到了很多数据(x, y)，那么什么样的模型才能让这些数据出现的可能性最大。这就是最大似然估计的出发点，也就是求解θ以最大化这些数据出现的概率，即最大化似然函数L(θ)。

关于最大似然估计方法更多解释可以看这里。

当然更多时候最大化的是logL(θ)，而不是直接最大化L(θ)，因为log函数单调递增函数，所以这个转化不会影响θ的最终取值。

l(θ)=logL(θ)=log∏i=1m12π??√σe?(y(i)?θTx(i))22σ2=∑i=1mlog12π??√σe?(y(i)?θTx(i))22σ2=mlog12π??√σ?1σ212∑i=1m(y(i)?θTx(i))2

因此最大化l(θ)也就是最小化：

12∑i=1m(y(i)?θTx(i))2

也就是之前出现的J(θ)。我们从概率和最大似然估计的角度解释了J(θ)选择这个二次式是合理的。