今天看啥 热点：
1.理论基础
& & 线性回归(Linear Regression)问题属于监督学习（Supervised Learning）范畴，又称分类（Classification）或归纳学习(Inductive Learning);这类分析中训练数据集中给出的数据类标是确定的；机器学习的目标是，对于给定的一个训练数据集，通过不断的分析和学习产生一个联系属性集合和类标集合的分类函数(Classification
Function)或预测函数(Prediction Function),这个函数称为分类模型(Classification Model)或预测模型(Prediction Model)；通过学习得到的模型可以是一个决策树,规格集,贝叶斯模型或一个超平面。通过这个模型可以对输入对象的特征向量预测或对对象的类标进行分类。
&& 回归问题中通常使用最小二乘(Least Squares)法来迭代最优的特征中每个属性的比重，通过损失函数(Loss Function)或错误函数(Error Function)定义来设置收敛状态，即作为剃度下降算法的逼近参数因子。
2.矩阵向量运算库jblas介绍
& & 由于spark MLlib中使用jlbas的线性代数运算库，因此学习和掌握jlbas库中基本的运算，对分析和学习spark中MLlib很多算法很有帮助；下面使用jlbas中DoubleMatrix矩阵对jlbas中基本运算进行简单介绍：
val matrix1 = DoubleMatrix.ones(10, 1) //创建所有值为1的10*1矩阵
val matrix2 = DoubleMatrix.zeros(10, 1) //创建所有值为0的10*1矩阵
matrix1.put(1, -10)
val absSum = matrix1.norm1() //绝对值之和
val euclideanNorm = matrix1.norm2() //欧几里德距离
val matrix3 = (matrix1.addi(matrix2))
val matrix4 = new DoubleMatrix(1, 10, (1 to 10).map(_.toDouble): _*) //创建Double向量对象
println(&print init value:matrix3=& + matrix3)
println(&print init value:matrix4=& + matrix4)
println(&matrix sub matrix:& + matrix3.sub(matrix4) + &,& + matrix4.sub(10)) //减法运算
println(&matrix add matrix:& + matrix3.add(matrix4) + &,& + matrix4.add(10)) //加法运算
println(&matrix mul matrix:& + matrix3.mul(matrix4) + &,& + matrix4.mul(10)) //乘法运算
println(&matrix div matrix:& + matrix3.div(matrix4) + &,& + matrix4.div(10)) //除法运算
println(&matrix dot matrix:& + matrix3.dot(matrix4)) //向量积
val matrix5 = DoubleMatrix.ones(10, 10)
println(&N*M Vector Matrix sub OP:\n& + matrix5.subRowVector(matrix4) + &\n& + matrix5.subColumnVector(matrix4)) //多对象减法运算
println(&N*M Vector Matrix add OP:\n& + matrix5.addRowVector(matrix4) + &\n& + matrix5.addColumnVector(matrix4)) //多对象加法运算
println(&N*M Vector Matrix mul OP:\n& + matrix5.mulRowVector(matrix4) + &\n& + matrix5.mulColumnVector(matrix4)) //多对象乘法运算
println(&N*M Vector Matrix div OP:\n& + matrix5.divRowVector(matrix4) + &\n& + matrix5.divColumnVector(matrix4)) //多对象除法运算
3.梯度下降(Gradient Descent)算法介绍
& & 梯度下降算法用于在迭代过程中逐渐降阶，不断更新特征权重向量，从而得到无限接近或拟合的最优特征权重向量 ；梯度下降算法主要有两种,第一种是批量梯度下降(Batch Gradient Descent)算法，此种方式实现过程是对权重向量进行累加，然后批量更新的一种方式，一般不实用于大规模数据集处理；另外一种是随机梯度下降(Stochastic
Gradient Descent)算法，这种方式对给定训练数据集中每个对象都进行权重计算和更新，在某些情况下容易收敛到局部最优解上。Spark MLlib库中主要使用随机梯度下降算法。为了更好的理解MLlib库中随机梯度算法(MLlib库中类后缀名以SGD结尾的所有算法)实现，下面是使用随机梯度算法进行直线拟合的Demo：
def sgdDemo {
val featuresMatrix: List[List[Double]] = List(List(1, 4), List(2, 5), List(5, 1), List(4, 2))//特征矩阵
val labelMatrix: List[Double] = List(19, 26, 19, 20)//真实值向量
var theta: List[Double] = List(0, 0)
var loss: Double = 10.0
i &- 0 until 1000 //迭代次数
if (loss & 0.01) //收敛条件loss&=0.01
var error_sum = 0.0 //总误差
var j = i % 4
var h = 0.0
for (k &- 0 until 2) {
h += featuresMatrix(j)(k) * theta(k)
} //计算给出的测试数据集中第j个对象的计算类标签
error_sum = labelMatrix(j) - h //计算给出的测试数据集中类标签与计算的类标签的误差值
var cacheTheta: List[Double] = List()
for (k &- 0 until 2) {
val updaterTheta = theta(k) + 0.001 * (error_sum) * featuresMatrix(j)(k)
cacheTheta = updaterTheta +: cacheTheta
} //更新权重向量
cacheTheta.foreach(t =& print(t + &,&))
print(error_sum + &\n&)
theta = cacheTheta
//更新误差率
var currentLoss: Double = 0
for (j &- 0 until 4) {
var sum = 0.0
for (k &- 0 until 2) {
sum += featuresMatrix(j)(k) * theta(k)
currentLoss += (sum - labelMatrix(j)) * (sum - labelMatrix(j))
loss = currentLoss
println(&loss-&&&&& + loss / 4 + &,i&&&&&& + i)
4.MLlib线性回归源码分析
& & &MLlib中可用的线性回归算法有：LinearRegressionWithSGD，RidgeRegressionWithSGD，LassoWithSGD；MLlib回归分析中涉及到的主要类有，GeneralizedLinearAlgorithm，GradientDescent。下面以对LinearRegressionWithSGD实现进行主要分析。
& & 第一步：在使用LinearRegressionWithSGD算法之前首先将输入数据解析成包含类标和特征向量的LabeledPoint对象的RDD弹性分布式数据集合。
& & 第二步：调用LinearRegressionWithSGD伴生对象中的train方法传输第一步创建的RDD集合和最大迭代次数，在train中主要实现创建一个新的LinearRegressionWithSGD对象，初始化梯度下降算使用使用最小二乘梯度下降算法SquaredGradient以及更新权重向量使用SimpleUpdater，执行父类GeneralizedLinearAlgorithm中的run方法进行权重向量和拦截参数计算，并返回训练得到的模型属性权重向量
& & LinearRegressionWithSGD伴生对象中train方法实现
def train(
input: RDD[LabeledPoint],
numIterations: Int,
stepSize: Double,//默认步长为1
miniBatchFraction: Double)//每次跌打使用的batch因子，默认为1
: LinearRegressionModel =
new LinearRegressionWithSGD(stepSize, numIterations, miniBatchFraction).run(input)
& &LinearRegressionWithSGD中run方法实现
def run(input: RDD[LabeledPoint], initialWeights: Array[Double]) : M = {
// Check the data properties before running the optimizer
if (validateData && !validators.forall(func =& func(input))) {//预验证输入数据的合法性，validators中存储验证的所有方法
throw new SparkException(&Input validation failed.&)
// Prepend an extra variable consisting of all 1.0's for the intercept.
val data = if (addIntercept) {//判读是否需要添加intercept
input.map(labeledPoint =& (labeledPoint.label, 1.0 +: labeledPoint.features))
input.map(labeledPoint =& (labeledPoint.label, labeledPoint.features))
}//对象转化为元组(类标签,特征)
val initialWeightsWithIntercept = if (addIntercept) {
0.0 +: initialWeights
initialWeights
}//初始化权重特征
val weightsWithIntercept = optimizer.optimize(data, initialWeightsWithIntercept)//返回优化后的权重
val (intercept, weights) = if (addIntercept) {
(weightsWithIntercept(0), weightsWithIntercept.tail)
(0.0, weightsWithIntercept)
logInfo(&Final weights & + weights.mkString(&,&))
logInfo(&Final intercept & + intercept)
createModel(weights, intercept)//使用计算后的权重向量和截距创建模型
其中optimizer.optimize(data, initialWeightsWithIntercept)是LinearRegressionWithSGD实现的核心，oprimizer的类型为GradientDescent，optimize方法中主要调用GradientDescent伴生对象的runMiniBatchSGD方法，返回当前迭代产生的最优特征权重向量。
GradientDescentd对象中optimize方法实现
def optimize(data: RDD[(Double, Array[Double])], initialWeights: Array[Double])
: Array[Double] = {
//返回优化后的权重向量，和迭代过程中误差
val (weights, stochasticLossHistory) = GradientDescent.runMiniBatchSGD(
numIterations,
miniBatchFraction,
initialWeights)
GradientDescent伴生对象中runMiniBatchSGD方法实现
def runMiniBatchSGD(
data: RDD[(Double, Array[Double])],
gradient: Gradient,//SquaredGradient—平方剃度下降算法
updater: Updater,//SimpleUpdater
stepSize: Double,//1.0
numIterations: Int,//100
regParam: Double,//0.0
miniBatchFraction: Double,//1.0
initialWeights: Array[Double]) : (Array[Double], Array[Double]) = {
val stochasticLossHistory = new ArrayBuffer[Double](numIterations)
val nexamples: Long = data.count()
val miniBatchSize = nexamples * miniBatchFraction
// Initialize weights as a column vector／／创建一维向量，第一个参数为行数，第二个参数为列数，第三个参数开始为值
var weights = new DoubleMatrix(initialWeights.length, 1, initialWeights:_*)
var regVal = 0.0
for (i &- 1 to numIterations) {
* 使用平方梯度下降算法
* gradientSum:本次选择迭代样本的梯度总和，
* lossSum:本次选择迭代样本的误差总和
val (gradientSum, lossSum) = data.sample(false, miniBatchFraction, 42+i).map {
case (y, features) =&//(标签，特征)
val featuresCol = new DoubleMatrix(features.length, 1, features:_*)
val (grad, loss) = pute(featuresCol, y, weights)//(特征，标签，特征属性权重向量)
* class SquaredGradient extends Gradient {
override def compute(data: DoubleMatrix, label: Double, weights: DoubleMatrix):
(DoubleMatrix, Double) = {
val diff: Double = data.dot(weights) - label//计算当前计算对象的类标签与实际类标签值之差
val loss = 0.5 * diff * diff//当前平方梯度下降值
val gradient = data.mul(diff)
(gradient, loss)
(grad, loss)//当前训练对象的特征权重向量和误差
}.reduce((a, b) =& (a._1.addi(b._1), a._2 + b._2)) //计算本次迭代所选训练数据新的特征权重向量之和和误差总和
* NOTE(Xinghao): lossSum is computed using the weights from the previous iteration
* and regVal is the regularization value computed in the previous iteration as well.
stochasticLossHistory.append(lossSum / miniBatchSize + regVal)//miniBatchSize=样本中对象数量＊batch因子，regVal为回归因子
val update = pute(
weights, gradientSum.div(miniBatchSize), stepSize, i, regParam)//weights：属性向量中设置的权重因子，regParam：为回归参数，stepSize：计算步长，i：当前迭代次数
* class SimpleUpdater extends Updater {
* weihtsOld:上一次迭代计算后的特征权重向量
* gradient:本次迭代计算的特征权重向量
* stepSize:迭代步长
* iter:当前迭代次数
* regParam:回归参数
override def compute(weightsOld: DoubleMatrix, gradient: DoubleMatrix,
stepSize: Double, iter: Int, regParam: Double): (DoubleMatrix, Double) = {
val thisIterStepSize = stepSize / math.sqrt(iter)//以当前迭代次数的平方根的倒数作为本次迭代趋近(下降)的因子
val normGradient = gradient.mul(thisIterStepSize)
(weightsOld.sub(normGradient), 0)//返回本次剃度下降后更新的特征权重向量
weights = update._1
regVal = update._2//使用SimpleUpdater值为0
logInfo(&GradientDescent finished. Last 10 stochastic losses %s&.format(
stochasticLossHistory.takeRight(10).mkString(&, &)))
(weights.toArray, stochasticLossHistory.toArray)
& 在MiniBatchSGD中主要实现对输入数据集进行迭代抽样，通过使用SquaredGradient作为梯度下降算法，使用SimpleUpdater作为更新算法，不断对抽样数据集进行迭代计算从而找出最优的特征权重向量解。
& 使用官方测试代码如下：
def linearRegressionAPITest(sc: SparkContext) {
val url = &/Users/yangguo/hadoop/spark/mllib/data/ridge-data/lpsa.data&
val data = sc.textFile(url)
val parseData = data.map { line =&
val parts = line.split(',')
LabeledPoint(parts(0).toDouble, parts(1).split(' ').map(x =& x.toDouble).toArray)
val numIterations = 20
val model = LinearRegressionWithSGD.train(parseData, numIterations)
val valuesAndPreds = parseData.map { point =&
val prediction = model.predict(point.features)
(point, prediction)
valuesAndPreds.foreach { case (v, p) =& print(&[& + v.label + &,& + p + &]&); v.features.foreach(base =& print(base + &--&)); println(&\n&) }
val isSuccessed = valuesAndPreds.map { case (v, p) =& math.pow((p - v.label), 2) }.reduce(_ + _) / valuesAndPreds.count
println(isSuccessed)
[1].http://rdc.taobao.org/?p=2163
[2].http://cs229.stanford.edu/notes/cs229-notes1.pdf
[3]..cn/s/blog_5jyq.html
[4].http://blog.csdn.net/pennyliang/article/details/6998517
[5].http://en.wikipedia.org/wiki/Lasso_(statistics)#Lasso_method
暂无相关文章
相关搜索：
相关阅读：
相关频道：
&&&&&&&&&&&&&&&&
WEB编程教程最近更新&|&&|&&|&&|&&
当前位置： >
[Machine Learning] 梯度下降（BGD）、随机梯度下降（SGD）、Mini-batch Gradient Descent、带Mini-batch的SGD
作者：互联网 & 来源：转载 &
浏览次数：
摘要: 一、回归函数及目标函数
以均方误差作为目标函数（损失函数），目的是使其值最小化，用于优化上式。
二、优化方式（Gradient Descent） 1、最速梯度下降法
也叫批量梯度下降法Batch Gradient Descent，BSD
a、对目标函数求导
b、沿导数相反方向移动theta
原因： （1）对于目标函数，theta的移动量应当如下，其中a为步长，p为方向向量
一、回归函数及目标函数
以均方误差作为目标函数（损失函数），目的是使其值最小化，用于优化上式。
二、优化方式（Gradient Descent）
1、最速梯度下降法
也叫批量梯度下降法Batch Gradient Descent，BSD
a、对目标函数求导
b、沿导数相反方向移动theta
（1）对于目标函数，theta的移动量应当如下，其中a为步长，p为方向向量。
（2）对J（theta）做一阶泰勒级数展开：
（3）上式中，ak是步长，为正数，可知要使得目标函数变小，则应当&0，并且其绝对值应当越大越好，这样下降的速度更快。在泰勒级数中，g代表J（theta k）的梯度，所以为了使得为负并且绝对值最大，应当使theta的移动方向与梯度g相反。
2、随机梯度下降法（stochastic gradient descent，SGD）
SGD是最速梯度下降法的变种。
使用最速梯度下降法，将进行N次迭代，直到目标函数收敛，或者到达某个既定的收敛界限。每次迭代都将对m个样本进行计算，计算量大。
为了简便计算，SGD每次迭代仅对一个样本计算梯度，直到收敛。伪代码如下（以下仅为一个loop，实际上可以有多个这样的loop，直到收敛）：
（1）由于SGD每次迭代只使用一个训练样本，因此这种方法也可用作online learning。
（2）每次只使用一个样本迭代，若遇上噪声则容易陷入局部最优解。
3、Mini-batch Gradient Descent
（1）这是介于BSD和SGD之间的一种优化算法。每次选取一定量的训练样本进行迭代。
（2）从公式上似乎可以得出以下分析：速度比BSD快，比SGD慢；精度比BSD低，比SGD高。
4、带Mini-batch的SGD
（1）选择n个训练样本（n&m，m为总训练集样本数）
（2）在这n个样本中进行n次迭代，每次使用1个样本
（3）对n次迭代得出的n个gradient进行加权平均再并求和，作为这一次mini-batch下降梯度
（4）不断在训练集中重复以上步骤，直到收敛。
版权所有 IT知识库 CopyRight (C)
, All Rights Reserved.Deep learning中的优化方法
　　三种常见优化算法：SGD（随机梯度下降），LBFGS（受限的BFGS），CG（共轭梯度法）。
& & &1.SGD（随机梯度下降）
& & & &随机梯度下降(Stochastic Gradient Descent, SGD)是随机和优化相结合的产物，是一种很神奇的优化方法，属于梯度下降的一种，适用于大规模问题。
　　要想扯清楚它，还得先谈谈梯度下降。众所周知，每个优化问题都会有一个目标函数F(w)F(w)，梯度下降采用迭代的策略，从初始点w0w0开始，每次沿着目标函数在当前点的负梯度方向前进一段距离，即
wt+1=wt?ηt?F(wt)wt+1=wt?ηt?F(wt)
只要步长ηtηt设置合理，就可以得到一个单调递减的序列f(w0),f(w1),?f(w0),f(w1),?，直至不再下降即可得到最优解w?w?。
对于一般的优化问题，梯度下降可以找到局部最优解，对于凸优化问题，梯度下降可以得到全局最优解，下面我们只考虑凸优化问题。
　　考虑如下的目标函数
F(w)=Ei∼Dfi(w)F(w)=Ei∼Dfi(w)
其中每个fifi都是关于ww的凸函数，下标ii服从分布DD。由期望的线性性有
?F(w)=Ei∼D?fi(w)?F(w)=Ei∼D?fi(w)
显然当DD是取值很多的离散分布或是连续分布时，?F(w)?F(w)计算开销很大甚至根本无法计算，这个方法也就行不通了。但这样的问题在机器学习领域又很常见，比如感知机、SVM、LASSO的优化目标都可以写成如下的形式
λΩ(w)+1N∑i=1NL(w,xi)λΩ(w)+1N∑i=1NL(w,xi)
其中Ω(w)Ω(w)是关于ww的凸正则化项，L(w,xi)L(w,xi)是在样本xixi上的损失，这里DD就是在NN个离散点上概率均为1/N1/N的离散分布。因此为了克服梯度下降的这个弱势，随机梯度下降应运而生。
　随机梯度下降的想法很简单，就是不直接计算梯度的精确值，而是用梯度的无偏估计
?fi(w)?fi(w)代替之作为下降方向，即在t+1t+1轮随机挑选出下标ii作如下更新：
wt+1=wt?ηt?fi(wt)wt+1=wt?ηt?fi(wt)
那么肯定有人要问，这么简单靠谱么？可以证明在一定条件下，这样得到的序列f(w0),f(w1),?f(w0),f(w1),?中的最小值依期望收敛到f(w?)f(w?)。具体来说，设
ηt≥0,&∑t=0∞η2t&∞,&∑t=0∞ηt=∞ηt≥0,&∑t=0∞ηt2&∞,&∑t=0∞ηt=∞
并假设存在常数GG满足对于任意t,it,i有E[∥?fi(wt)∥2]≤G2E[∥?fi(wt)∥2]≤G2及常数RR满足E[∥w0?w?∥2]≤R2E[∥w0?w?∥2]≤R2，并记fbest(t)=min{f(w0),?,f(wt)}fbest(t)=min{f(w0),?,f(wt)}，那么当t→∞t→∞时，E[fbest(t)]→f(w?)E[fbest(t)]→f(w?)。
　　设t+1t+1轮随机挑出的下标为ii，那么
∥wt+1?w?∥2=∥wt?ηt?fi(w)?w?∥2=∥wt?w?∥2?2ηt?fi(w)?(wt?w?)+η2t∥?fi(w)∥2∥wt+1?w?∥2=∥wt?ηt?fi(w)?w?∥2=∥wt?w?∥2?2ηt?fi(w)?(wt?w?)+ηt2∥?fi(w)∥2
结合条件期望的线性性有
E[∥wt+1?w?∥2|wt]=E[∥wt?w?∥2|wt]?2ηtE[?fi(w)?(wt?w?)|wt]+η2tE[∥?fi(w)∥2|wt]=∥wt?w?∥2?2ηt?F(wt)(wt?w?)+η2tE[∥?fi(w)∥2|wt]≤∥wt?w?∥2?2ηt(F(wt)?F(w?))+η2tG2E[∥wt+1?w?∥2|wt]=E[∥wt?w?∥2|wt]?2ηtE[?fi(w)?(wt?w?)|wt]+ηt2E[∥?fi(w)∥2|wt]=∥wt?w?∥2?2ηt?F(wt)(wt?w?)+ηt2E[∥?fi(w)∥2|wt]≤∥wt?w?∥2?2ηt(F(wt)?F(w?))+ηt2G2
两边同时对wtwt取期望，由重期望公式
E[∥wt+1?w?∥2]≤E[∥wt?w?∥2]?2ηt(E[F(wt)]?F(w?))+η2tG2E[∥wt+1?w?∥2]≤E[∥wt?w?∥2]?2ηt(E[F(wt)]?F(w?))+ηt2G2
重复利用该式可得
E[∥wt+1?w?∥2]≤E[∥w0?w?∥2]?2∑j=0tηj(E[F(wj)]?F(w?))+G2∑j=0tη2jE[∥wt+1?w?∥2]≤E[∥w0?w?∥2]?2∑j=0tηj(E[F(wj)]?F(w?))+G2∑j=0tηj2
注意E[∥wt+1?w?∥2]≥0E[∥wt+1?w?∥2]≥0以及E[∥w0?w?∥2]≤R2E[∥w0?w?∥2]≤R2，于是
2∑j=1tηj(E[F(wj)]?F(w?))≤R2+G2∑j=0tη2j2∑j=1tηj(E[F(wj)]?F(w?))≤R2+G2∑j=0tηj2
结合E[Fbest(t)]≤E[F(wj)]E[Fbest(t)]≤E[F(wj)]可知
E[Fbest(t)]?F(w?)≤R2+G2∑tj=0η2j2∑tj=1ηjE[Fbest(t)]?F(w?)≤R2+G2∑j=0tηj22∑j=1tηj
由于∑∞t=1ηt=∞∑t=1∞ηt=∞，故当t→∞t→∞时有E[Fbest(t)]→F(w?)E[Fbest(t)]→F(w?)。
　　此外，由Markov不等式知对于??&0??&0有
P(Fbest(t)?F(w?)≥?)≤E[Fbest(t)?F(w?)]?≤R2+G2∑tj=0η2j2?∑tj=1ηjP(Fbest(t)?F(w?)≥?)≤E[Fbest(t)?F(w?)]?≤R2+G2∑j=0tηj22?∑j=1tηj
即当t→∞t→∞时有P(Fbest(t)?F(w?)≥?)→0P(Fbest(t)?F(w?)≥?)→0。
下面举几个机器学习里的例子，设训练集有NN个样本{(x1,y1),(x2,y2),…,(xN,yN)}{(x1,y1),(x2,y2),…,(xN,yN)}，于是此时DD就是在NN个离散点上概率均为1/N1/N的离散分布，易知有
?F(w)=?(Ω(w)+1N∑i=1NL(w,xi))=1N∑i=1N?(Ω(w)+L(w,xi))=E[?(Ω(w)+L(w,xi))]?F(w)=?(Ω(w)+1N∑i=1NL(w,xi))=1N∑i=1N?(Ω(w)+L(w,xi))=E[?(Ω(w)+L(w,xi))]
于是随机梯度下降就是每次随机选取一个样本xixi，以??(Ω(w)+L(w,xi))??(Ω(w)+L(w,xi))作为下降方向。
感知机可形式化成如下的优化问题
minw&&&1N∑i=1Nmax{0,?yiw?xi}minw&&&1N∑i=1Nmax{0,?yiw?xi}
设t+1t+1轮随机挑出的样本为(xi,yi)(xi,yi)，那么对应的更新公式为
wt+1=wt+ηt{yixi0yiw?txi&0otherwisewt+1=wt+ηt{yixiyiwt?xi&00otherwise
SVM可形式化成如下的优化问题
minw&&&λ2∥w∥2+1N∑i=1Nmax{0,1?yiw?xi}minw&&&λ2∥w∥2+1N∑i=1Nmax{0,1?yiw?xi}
设t+1t+1轮随机挑出的样本为(xi,yi)(xi,yi)，那么对应的更新公式为
wt+1=wt?ηt{λwt?yixiλwtyiw?txi&1otherwisewt+1=wt?ηt{λwt?yixiyiwt?xi&1λwtotherwise
LASSO可形式化成如下的优化问题
minw&&&λ∥w∥1+1N∑i=1N12(w?xi?yi)2minw&&&λ∥w∥1+1N∑i=1N12(w?xi?yi)2
设w=u?vw=u?v且u≥0,v≥0u≥0,v≥0，ee为全11向量，优化问题可重写为
minu,v&&&λu?e+λv?e+1n∑i=1N12(u?xi?v?xi?yi)2minu,v&&&λu?e+λv?e+1n∑i=1N12(u?xi?v?xi?yi)2
设t+1t+1轮随机挑出的样本为(xi,yi)(xi,yi)，那么对应的更新公式为
ut+1=max{0,ut?ηt(λe+(w?txi?yi)xi)}vt+1=max{0,vt?ηt(λe?(w?txi?yi)xi)}ut+1=max{0,ut?ηt(λe+(wt?xi?yi)xi)}vt+1=max{0,vt?ηt(λe?(wt?xi?yi)xi)}
　　最后再提一个小，以支持向量机为例，它的更新公式为
wt+1=(1?ληt)wt+ηt{yixi0yiw?txi&1otherwise(1)(1)wt+1=(1?ληt)wt+ηt{yixiyiwt?xi&10otherwise
当xx维度很高而非零很少时，+yixi+yixi可以很高效地算出来，但是第一项(1?ληt)wt(1?ληt)wt的计算代价就有点高了，因为ww一般来说不是稀疏的，一个就是做个变量代换
wt=utatwt=utat
其中αtαt是标量，于是式(11)可以转化为如下只涉及标量计算和稀疏向量操作的更新过程
ztat+1ut+1=yiu?txi/at=at1?ληt=ut+at+1ηt{yixi0zt&1otherwise
　　SGD优点：实现简单，当训练样本足够多时优化速度非常快。
　　SGD缺点：需要人为调整很多参数，比如学习率，收敛准则等。另外，它是序列的方法，不利于GPU并行或分布式处理。
2.L-BFGS（受限的BFGS）
&& L-BFGS即Limited-memory BFGS，在之前的BFGS算法中，我们可以不矩阵，而是存储最近次迭代
&& 的曲率信息，即和。当完成一次迭代后，最旧的一次曲率的信息将被删除，而最新的曲率将被保存下来，所
&& 以这样就保证了保存的曲率信息始终都来自最近的次迭代。在实际工程中取3到20之间的值效果比较好。
&& 在之前的BFGS算法中，有如下公式
&& 那么同样有
&& 将这个式子带入，得到
&&&整理一下，一直递推次下去，就有
&&&每次迭代初始值的设定，在实践中常用的方法是
&& 利用最近一次的曲率信息来估计真实Hessian矩阵的大小，这样使得当前搜索方向较为理想，不至于跑得太偏。
3.CG（共轭梯度法）
共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一.
& &最初是由计算数学家Hestenes和几何学家Stiefel于1952年为求正定系数矩阵线性方程组而独立提出的.他们合作的著名Method of conjugate gradients for solving linear s 被认为是共轭梯度法的奠基性文章。1964年，Fletcher和Reeves将此方法推广到非线性最优化，
得到了求解一般函数极小值的共轭梯度法.共轭梯度法的收敛性分析的早期工作主要由Fletcher、 Powell、Beale等学者给出.
(1) 建立在二次模型上，具有二次终止性．
(2) 一种有效的算法，克服了最速下降法的锯齿现象，又避免了牛顿法的计算量大和局部收敛性的缺点．
(3) 算法简单，易于，无需计算二阶导数，存储 小等优点，是求解中等规模优化问题的主要方法．
& &共轭梯度法（conjugate gradient method, CG）是以共轭方向（conjugate direction）作为搜索方向的一类算法。CG法是由Hesteness和Stiefel于1952年为求解线性方程组而提出的。后来用于求解无约束最优化问题，它是一种重要的数学优化方法。这种方法具有二次终止性。
& &CG的基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿着此组方向进行搜索，求出目标函数的极小点。
　　各种deep learning中常见方法（比如说Autoencoder，RBM，DBN，ICA，Sparse coding）的区别是：目标函数形式不同。这其实才是最本质的区别，由于目标函数的不同导致了对其优化的方法也可能会不同，比如说RBM中目标函数跟能量有关，采用CD优化的，而Autoencoder目标函数为理论输出和实际输出的MSE，由于此时的目标函数的偏导可以直接被计算，所以可以用LBFGS，CG等方法优化，其它的类似。所以不能单从网络的结构来判断其属于Deep learning中的哪种方法，比如说我单独给定64-100的2层网络，你就无法知道它属于deep
learning中的哪一种方法，因为这个网络既可以用RBM也可以用Autoencoder来训练。
　　通过实验得出的结论是：不同的优化算法有不同的优缺点，适合不同的场合，比如LBFGS算法在参数的维度比较低（一般指小于10000维）时的效果要比SGD（随机梯度下降）和CG（共轭梯度下降）效果好，特别是带有convolution的模型。而针对高维的参数问题，CG的效果要比另2种好。也就是说一般情况下，SGD的效果要差一些，这种情况在使用GPU加速时情况一样，即在GPU上使用LBFGS和CG时，优化速度明显加快，而SGD算法优化速度提高很小。在单核上，LBFGS的优势主要是利用参数之间的2阶近视特性来加速优化，而CG则得得益于参数之间的共轭信息，需要计算器Hessian矩阵。
　　不过当使用一个大的minibatch且采用线搜索的话，SGD的优化性能也会提高。
　　在单核上比较SGD，LBFGS，CG三种算法的优化性能，当针对Autoencoder模型。结果如下：
　　可以看出，SGD效果最差。
　　同样的情况下，训练的是Sparse autoencoder模型的比较情况如下：
　　这时SGD的效果更差。这主要原因是LBFGS和CG能够使用大的minibatch数据来估算每个节点的期望激发值，这个值是可以用来约束该节点的稀疏特性的，而SGD需要去估计噪声信息。
　　当然了作者还有在GUP，convolution上也做了不少实验。
　　最后，作者训练了一个2隐含层（这2层不算pooling层）的Sparse autocoder网络，并应用于MNIST上，其识别率结果如下：
下面是作者code主要部分的一些注释：
optimizeAutoencoderLBFGS.m(实现deep autoencoder网络的参数优化过程):
function [] = optimizeAutoencoderLBFGS(layersizes, datasetpath, ...
finalObjective)
% train a deep autoencoder with variable hidden sizes
% layersizes : the sizes of the hidden layers. For istance, specifying layersizes =
[200 100] will create a
looks like input -& 200 -& 100 -& 200
-& output (same size as input). Notice the mirroring structure of the
autoencoders. Default layersizes = [2*3072 100]
% datasetpath: the path to the CIFAR dataset (where we find the *.mat
files). see loadData.m
% finalObjective: the final objective that you use to compare to
terminate your optimization. To qualify, the objective
function on the entire training set must be below this
% Author: Quoc V. Le (quocle@stanford.edu)
%% Handle default parameters
if nargin & 3 || isempty(finalObjective)
finalObjective = 70; % i am just making this up, the evaluation objective
% will be much lower
if nargin & 2 || isempty(datasetpath)
datasetpath = '.';
if nargin & 1 || isempty(layersizes)
layersizes = [2*3072 100];
layersizes = [200 100];
%% Load data
loadData %traindata 3072*10000的，每一列表示一个向量
%% Random initialization
initializeW%看作者对应该部分的code，也没有感觉出convolution和pooling的影响啊，怎么它就连接起来了呢
%% Optimization: minibatch L-BFGS
% Q.V. Le, J. Ngiam, A. Coates, A. Lahiri, B. Prochnow, A.Y. Ng.
% On optimization methods for deep learning. ICML, 2011
addpath minFunc/
options.Method = 'lbfgs';
options.maxIter = 20;
options.display = 'on';
options.TolX = 1e-3;
perm = randperm(size(traindata,2));
traindata = traindata(:,perm);% 将训练样本随机排列
batchSize = 1000;%因为总共样本数为10000个，所以分成了10个批次
maxIter = 20;
for i=1:maxIter
startIndex = mod((i-1) * batchSize, size(traindata,2)) + 1;
fprintf('startIndex = %d, endIndex = %d\n', startIndex, startIndex + batchSize-1);
data = traindata(:, startIndex:startIndex + batchSize-1);
[theta, obj] = minFunc( @deepAutoencoder, theta, options, layersizes, ...
if obj &= finalObjective % use the minibatch obj as a heuristic for stopping
% because checking the entire dataset is very
% expensive
% yes, we should check the objective for the entire training set
trainError = deepAutoencoder(theta, layersizes, traindata);
if trainError &= finalObjective
% now your submission is qualified
%% write to text files so that we can test your program
writeToTextF
deepAutoencoder.m:（深度网络代价函数及其导数的求解函数）:
function [cost,grad] = deepAutoencoder(theta, layersizes, data)
% cost and gradient of a deep autoencoder
% layersizes is a vector of sizes of hidden layers, e.g.,
% layersizes[2] is the size of layer 2
% this does not count the visible layer
% data is the input data, each column is an example
% the activation function of the last layer is linear, the activation
% function of intermediate layers is the hyperbolic tangent function
% WARNING: the code is optimized for ease of implemtation and
% understanding, not speed nor space
%% FORCING THETA TO BE IN MATRIX FORMAT FOR EASE OF UNDERSTANDING
% Note that this is not optimized for space, one can just retrieve W and b
% on the fly during forward prop and backprop. But i do it here so that the
% readers can understand what's going on
layersizes = [size(data,1) layersizes];
l = length(layersizes);
for i=1:l-1
lold = lnew + 1;
lnew = lnew + layersizes(i) * layersizes(i+1);
W{i} = reshape(theta(lold:lnew), layersizes(i+1), layersizes(i));
lold = lnew + 1;
lnew = lnew + layersizes(i+1);
b{i} = theta(lold:lnew);
% handle tied-weight stuff
for i=l:2*(l-1)
lold = lnew + 1;
lnew = lnew + layersizes(l-j);
W{i} = W{l - j}'; %直接用encoder中对应的转置即可
b{i} = theta(lold:lnew);
j = j + 1;
assert(lnew == length(theta), 'Error: dimensions of theta and layersizes do not match\n')
%% FORWARD PROP
for i=1:2*(l-1)-1
[h{i} dh{i}] = tanhAct(bsxfun(@plus, W{i}*data, b{i}));
[h{i} dh{i}] = tanhAct(bsxfun(@plus, W{i}*h{i-1}, b{i}));
h{i+1} = linearAct(bsxfun(@plus, W{i+1}*h{i}, b{i+1}));
%% COMPUTE COST
diff = h{i+1} -
M = size(data,2);
cost = 1/M * 0.5 * sum(diff(:).^2);% 纯粹标准的autoencoder，不加其它比如sparse限制
%% BACKPROP
if nargout & 1
outderv = 1/M *
for i=2*(l-1):-1:2
Wgrad{i} = outderv * h{i-1}';
bgrad{i} = sum(outderv,2);
outderv = (W{i}' * outderv) .* dh{i-1};
Wgrad{1} = outderv * data';
bgrad{1} = sum(outderv,2);
% handle tied-weight stuff
for i=l:2*(l-1)
Wgrad{l-j} = Wgrad{l-j} + Wgrad{i}';
j = j + 1;
% dump the results to the grad vector
grad = zeros(size(theta));
for i=1:l-1
lold = lnew + 1;
lnew = lnew + layersizes(i) * layersizes(i+1);
grad(lold:lnew) = Wgrad{i}(:);
lold = lnew + 1;
lnew = lnew + layersizes(i+1);
grad(lold:lnew) = bgrad{i}(:);
for i=l:2*(l-1)
lold = lnew + 1;
lnew = lnew + layersizes(l-j);
grad(lold:lnew) = bgrad{i}(:);
j = j + 1;
%% USEFUL ACTIVATION FUNCTIONS
function [a da] = sigmoidAct(x)
a = 1 ./ (1 + exp(-x));
if nargout & 1
da = a .* (1-a);
function [a da] = tanhAct(x)
a = tanh(x);
if nargout & 1
da = (1-a) .* (1+a);
function [a da] = linearAct(x)
if nargout & 1
da = ones(size(a));
initializeWeights.m（参数初始化赋值，虽然是随机，但是有一定要求）:
%% Random initialization
% X. Glorot, Y. Bengio.
% Understanding the dif铿乧ulty of training deep feedforward neural networks.
% AISTATS 2010.
% QVL: this initialization method appears to perform better than
% theta = randn(d,1);
s0 = size(traindata,1);% s0涓烘牱鏈殑缁存暟
layersizes = [s0 layersizes];%输入层-hidden1-hidden2，这里是3072-6144-100
l = length(layersizes);%缃戠粶涓殑灞傛暟锛屼笉鍖呭惈瑙ｇ爜閮ㄥ垎锛屽鏋滄槸2涓殣鍚眰鐨勮瘽锛岃繖閲宭=3
for i=1:l-1%1到3之间
lold = lnew + 1;
lnew = lnew + layersizes(i) * layersizes(i+1);
= sqrt(6) / sqrt(layersizes(i+1)+layersizes(i));
A = rand(layersizes(i+1), layersizes(i))*2*r - %reshape(theta(lold:lnew), layersizes(i+1), layersizes(i));
theta(lold:lnew) = A(:); %相当于权值W的赋值
lold = lnew + 1;
lnew = lnew + layersizes(i+1);
A = zeros(layersizes(i+1),1);
theta(lold:lnew) = A(:);%相当于偏置值b的赋值
end %以上是encoder部分
for i=l:2*(l-1) %1到4之间，下面开始decoder部分
lold = lnew + 1;
lnew = lnew + layersizes(l-j);
theta(lold:lnew)= zeros(layersizes(l-j),1);
j = j + 1;
theta = theta';
layersizes = layersizes(2:end); %去除输入层