梯度下降(Gradient Descent)，一句代码，一个式子

文章来源：http://www.zhizhihu.com/html/y2011/3632.html

一直以来，总是觉得国外的PhD们的教育以及课程的安排很好很强大，虽然是说很累作业多工作量大，但是功率大了，效果好点儿，浪费的时间也少，年轻人哪有怕苦怕累的。比比身边好多每天睡超过12小时的研究生们，不知道是谁更幸福一点儿。

我也经常拿我所在的大学的研究生博士跟自己所了解的美国那边的phd比比，说实话感觉总是有点硬性的差距在里面。我所处在大学的研究生教育，就拿安排的课程来说，觉得真是有用的不多，虽然也有一些好课，但是鉴于教学质量的问题，往往又起不到好的效果，最重要的是，这里的学生真是太轻松安逸了，不会有有用的作业，也不会有作业让你熬夜计算或者coding，也不会有作业让你在图书馆过夜，图书馆，哎，也是在按时上下班。

国外大牛的课程设计的也不错，拿美女教授Kristen Grauman举个例子，你可以去看看她得研究生以及本科的课程，看着都羡慕啊，是吧。Andrew
Ng的课也不用说了，还有很多很多自身修养不错的老师的课设计的都不错。及时你对这个方向开始不是很懂，那么这一门课下来，读了这么多论文以及写了这么多作业和coding，加上一门课所辐射的一个领域，相信应该会学到很多东西，况且老师教学又那么认真，及时你懒，那么也会逼着你去学到了。

课程设计真是一种艺术，也是一种优化。学到的不仅仅是理论，也是一种理解和实现理论的方法以及过程的体验。

--------------------------------



我经常会有一种感觉，就是感觉在自己实现一些算法的时候总是觉得不顺利或者甚至不想去实现，一方面开源的工具很多，另一方面懒吧。

其实每个初学算法的人都应该去亲自体验一下才好，我觉得这个过程当中最重要的就是建立起书面论文理论到具体coding实现的一个映射，理论-实现，这个映射相信是很多人的瓶颈，要努力克服这个。

就拿梯度下降(Gradient Descent)这个小算法来举个例子吧，算法简单，对吧？但是忽然让你马上实现一下，能够5分钟迅速搞定么？3分钟？？允许用matlab的话，方便点儿是吧。

------------------------------------------------

即使以前看过梯度下降算法，但是真的要去用去实现的时候，还真是得重新再看看。针对具体问题吧，先选择个问题，回归问题。再找点儿数据，索性就把libsvm包里的例子数据heart_scale.mat拿来吧，可以自己去下载。虽然是分类的，但是也可以看做回归呗，有些原理也相似。数据中heart_scale_inst包括270个13维的样本，label全部是+-1，这里看做回归吧。

假设要学习这么一个函数：

h(x)=hθ(x)=θ0+θ1x1+θ2x2+⋯
那么损失函数可以定义成：

J(Θ)=12∥ΘX−Y∥2
那么我们的目标很简单，就是求损失的最小值的时候的解：

Θ∗=argminΘJ(Θ)
像这种优化问题有很多方法，那咱们先直接求导吧，对于求导过程，好多还是不理解，可以用这种方法：

首先定义损失变量：

ri=∑j=1nXijθj−yi
那么损失函数就可以表示成：

J=12∑i=1mr2i
一步一步的求导，先求∂ri∂θj的导数：

∂J∂θj=∑i=1mri∂ri∂θj(j=1,2,…,n)
再求：

∂ri∂θj=Xij
那么把分步骤合起来就是：

∂J∂θj=∑i=1m(∑k=1nXikθk−yi)Xij(j=1,2,…,n)
导数为0，那么有：

∑i=1m(∑k=1nXikθˆk−yi)Xij=0(j=1,2,…,n)
整理一下：

∑i=1m∑k=1nXijXikθˆk=∑i=1mXijyi(j=1,2,…,n)
用矩阵符号将上面的细节运算抽象一下：

∂J(Θ)∂Θ=XTXΘ−XTY=0
让导数为0，那么求得的解为：

Θ=(XTX)−1XTY
但是我们知道求矩阵(XTX)的逆复杂度有点儿高，O(n3)，如果n很大，那么也受不了。

可以用最小二乘或者梯度下降来求解，这里我们看看梯度下降的实现，梯度下降的思想不难，只要确定好梯度以及梯度的方向就ok，因为是梯度的反方向去下降，所以在对参数更新的时候要注意：

Θi=Θi−1−γ∇J(Θ)=Θi−1−γ∂J(Θ)∂Θ
其中γ就是下降的速度了，这个很敏感的，一般是一个小的数值，可以从0.01开始尝试，越大下降越快，收敛越快。当然下降的速率可以改成自适用的，就是根据梯度的强弱适当调整步伐，这样效果还好一点儿。



上图就是迭代目标函数值的情况，迭代终止的条件有很多种（比如检测J(Θ)是否变化），这里取得：

∥∥Θi−Θi−1∥∥<ε

代码也很简单，如果用matlab的话，矩阵计算就容易很多：

clc;
clear
% load data
heart_scale = load('heart_scale');
X = heart_scale.heart_scale_inst;
Y = heart_scale.heart_scale_label;
 
epsilon = 0.0003;
gamma= 0.0001;
 
w_old=zeros(size(X,2),1);
k=1;
figure(1);
while 1
minJ_w(k) = 1/2 * (norm(X*w_old - Y))^2;
w_new = w_old - gamma*(X'*X*w_old - X'*Y);
fprintf('The %dth iteration, minJ_w = %f, \n',k,minJ_w(k));
 
if norm(w_new-w_old) < epsilon
W_best = w_new;
break;
end
w_old = w_new;
k=k+1;
end
 
plot(minJ_w);

运行的结果：

The 1th iteration, minJ_w = 135.000000,

The 2th iteration, minJ_w = 128.785026,

The 3th iteration, minJ_w = 123.210569,

The 4th iteration, minJ_w = 118.202908,

The 5th iteration, minJ_w = 113.697504,

The 6th iteration, minJ_w = 109.637791,

The 7th iteration, minJ_w = 105.974133,

The 8th iteration, minJ_w = 102.662919,

The 9th iteration, minJ_w = 99.665786,

The 10th iteration, minJ_w = 96.948948,

The 11th iteration, minJ_w = 94.482602,

The 12th iteration, minJ_w = 92.240436,

The 13th iteration, minJ_w = 90.199178,

The 14th iteration, minJ_w = 88.338227,

The 15th iteration, minJ_w = 86.639312,

....................