1、这是75%完成度的SGD代码
主要的程序在八月中旬基本完成了，后边忙别的事情，实在懒得再仔细写了，目前的代码是75%完成度的，已知bug和可以改进的地方有：
1、在多次进行下降后，参数矩阵可能变成NaN值。（即由于参数计算导致除零错，目前尚未排查，估计在6%错误率左右的时候出现）
2、并发性可以进行改进。包括采用GPU进行计算的可能。
3、断点保存可以改进。目前保存每次的参数矩阵，可以从上次参数矩阵进行下降。但尚未保存每次退出时进行下降的训练例子位置等信息。
这个SGD的效果还行，应该不比原例子差，不过性能就差了不止十万八千里（囧）。但是考虑到这是第一次写这种算法，没有必要对自己吹毛求疵，而且目的是为了学习和理解SGD算法，自我安慰安慰。
原例子在http://deeplearning.net/tutorial/gettingstarted.html，采用python写的，主要数学工具是Theano。
这个SGD代码的程序一会儿我上载到资源区。
2、为什么要用 C++实现SGD算法
这是因为没有搞定Theano。Logistic-SGD所需环境包括Theano，但是在尝试在Win7/IronPython VS2010环境下安装Theano未果之后，只能采取手动挡来研究一下Logistic-SGD的代码了。
基本思路是通过IronPython+VS2010 .Net环境，在找一些数学库（用于矩阵运算等等），替换Theano的部分内容，使得程序能够跑起来。
结果后来VS2010.Net操作IronPython没有搞定numpy库（抓狂，在外边直接可以用numpy，在C++初始化IronPython却不能使用numpy），加上学的深入一点后，似乎不用特别复杂的数学库，直接实现也可以。于是就有了这个程序。
3、SGD理论
4、程序
原例子中的mnist.pkg.gz是压缩后的数据族，为了能够便于C++直接处理，我们首先把mnist.pkg.gz改成易于处理的数据-结果对，即：将mnist.pkg解压为6个文件：
train_set_x
train_set_y
valid_set_x
valid_set_y
test_set_x
test_set_y
对应train,valid,test三个集合（每个集合两对，即数据-结果对）。train_set_x中的784行对应train_set_y中的一行，即代表一个测试数据集合。其中train_set_x每一行是一个浮点数，748行（784为28x28的灰度浮点数）代表一个数字的灰度图，对应train_set_y的一行，即对应的数字。
这个python程序如下，很简单，写成文本文件好了，然后再写一个C++程序读这个文本文件，转成相关的 二进制文件。
import cPickle
f = open(‘mnist.pkg’, 'rb')
train_set,valid_set, test_set = cPickle.load(f)
f.close()
test_set_x, test_set_y = test_set
train_set_x, train_set_y = train_set
valid_set_x, valid_set_y = valid_set
fileHandle=open("test_set_x","w")
cnt = 0
for i in test_set_x:
ifcnt% 10 ==0:
printcnt
cnt= cnt +1
forj in i:
fileHandle.write("%.6f\n"%j)
fileHandle.close()
cnt = 0
fileHandle=open("test_set_y","w")
for i in test_set_y:
ifcnt% 10 ==0:
printcnt
cnt= cnt +1
fileHandle.write("%d\n"%i)
fileHandle=open("train_set_x","w")
cnt=0
for i in train_set_x:
ifcnt% 10 ==0:
printcnt
cnt= cnt +1
forj in i:
fileHandle.write("%.6f\n"%j)
fileHandle.close()
fileHandle=open("train_set_y","w")
cnt=0
for i in train_set_y:
ifcnt% 10 ==0:
printcnt
cnt= cnt +1
fileHandle.write("%d\n"%i)
为了方便，已经做了二进制文件。
4.1、程序结构
整个程序由三部分组成：主程序、SGD主程序、文件操作支持函数组成。
数据格式在dataset.h中定义，一个集合（训练、验证、测试）为：
|
/* DSET_X 为每一个样例的灰度图，由784个double构成 typedef double DSET_X[784];
typedef struct st_set_array { int m_iCount; DSET_X *pArrayX; /* 每一个样本的灰度图 */ int *pArrayY; }SET_ARRAY; |
训练对应的参数w矩阵和b向量定义如下：
|
/* * WMATRIX W matrix parameter * bVECTOR b vector parameter */ typedef DSET_X WMATRIX[10]; typedef double bVECTOR[10];
|
矩阵定义其实可以用简单的wMatrix[ ][ ]来描述，但是我比较笨，懒得去思考C++二维数组的行列，索性就简单的先定义DSET_X然后再定义一维数组，之后的操作中全部采用一维数组的操作，避免因为笨导致的问题。
定义了一个控制块，用于各个函数中传递参数使用，
|
typedef struct st_sgd { double m_dLearnRate; int m_iEpoch; int m_iBatchSize; /* ----------------------------- */ WMATRIX m_wMatrix; bVECTOR m_bVector; SET_ARRAY *m_pTrainSet; SET_ARRAY *m_pValidSet; SET_ARRAY *m_pTestSet; /* ----------------------------- */ MATRIXINDEX m_MatrixIndex; double *m_Pi; int m_iValidationFrequency; int m_iPatience; int m_iPatienceIncrease;
}ST_SGD;
|
这里的内容可以不求甚解。看变量名也应该可以猜到干啥的，m_dLearnRate是学习率（就是参数delta的系数），Epoch是最大学习迭代次数。
按照整体的思想，我每次训练一批（一个BatchSize个训练样本），就调整并dump一次参数，然后训练一定批次后，用validate样本进行验证，当valid有改善时，用test样本进行验证。
中间退出的话，由于保存了每批dump的参数，下次直接可以用最新的参数进行训练。
Ok，let‘s do it。
4.2、主程序
主程序很简单，在minst04.cpp中（为啥是04？因为这是我第四次重写这个程序）。关键点在以下几行代码：
1、首先是读取训练样本文件
LoadDataSet(&TrainSet,&TestSet,&ValidSet);
这个是从"C:\\DEVELOP\\MNIST\\“下读取相关的二进制训练集合（LoadDataSet的函数实现也在mnist04.cpp中，其具体读取文件的函数实现在FileConvert.cpp中。
2、初始化sgd控制结构
ZeroMemory(&m_sgd,sizeof(m_sgd));
m_sgd.m_dLearnRate=0.13; 学习率0.13
m_sgd.m_iBatchSize=600; 600个样本一组
m_sgd.m_iEpoch=1000; 最多训练1000次
m_sgd.m_iPatience=5000; 。。。难以言表
m_sgd.m_iPatienceIncrease=2; 。。。难以言表
m_sgd.m_pTestSet = &TestSet; 测试集合
m_sgd.m_pTrainSet = &TrainSet; 训练集合
m_sgd.m_pValidSet = &ValidSet; 验证集合
3、读取上次dump的参数集合。
我们每次dump的参数文件在”C:\\DEVELOP\\MNIST\\parms\\“目录下，文件名以parms打头，因此读取上次dump的参数是找到最新的parms文件，然后读取之：
GetLatestFile("C:\\DEVELOP\\MNIST\\parms\\","parms*.txt",sName);
LoadParameters(sName,m_sgd.m_wMatrix,m_sgd.m_bVector);
相关函数在FileConvert.cpp中。
读取完，马上dump一份（可以看看读的是否正确。。。-_-!）
DumpMatrix("C:\\DEVELOP\\MNIST\\parms\\",m_sgd.m_wMatrix,m_sgd.m_bVector);
4、开始SGD
test_sgd(&m_sgd);
4.3、SGD
SGD的核心就是算几个数：
第一个是算
，即，给定某一个样本X，在W、b参数下，得出结果是i的概率。这个简单，在Softmax.cpp中，我们实现了两个（其实是一个意思，一个是快速查表算法Pn，一个是标准算法PnNoFast），
以PnNoFast作为样例解剖，原型：
double PnNoFast(int iY,int iItem,SET_ARRAY *pSet,WMATRIX w,bVECTOR b)
返回pSet集合中的第iItem个样本，在w和b参数下，当输出分类是iY时的概率。
关键代码：
pw = &w[iY];
px = &pSet->pArrayX[iItem];
dExp1 = 0;
/* 计算矩阵相乘， Wi * X + bi */
for(i=0;i<784;i++)
{
dExp1 += (*pw)[i]*(*px)[i];
}
dExp1 += b[iY];
dUpperValue=exp(dExp1);
这段代码实现了
dLowerValue = 0;
for(j=0;j<10;j++)
{
dExp1 = 0;
pw = &w[j];
for(i=0;i<784;i++)
{
dExp1 +=(*pw)[i]*(*px)[i];
}
dExp1 += b[j];
/*sum e^(WjX+bj) */
dLowerValue += exp(dExp1);
}
这段代码实现了
当然dUpperValue/dLowerValue=
第二个是算偏导数。
double Dpi_dWmk(MATRIXINDEX pMatrixIndex,int iItem,int m,int k,WMATRIXwMatrix,bVECTOR bVector,SET_ARRAY *pSet)
{
/* dpi/dwmk = (u'v-uv')/v^2 */
double du,dv,u,v;
du =DuI(pMatrixIndex,iItem,m,k,wMatrix,bVector,pSet);
dv =Dvi(pMatrixIndex,iItem,m,k,wMatrix,bVector,pSet);
u =Ui(pMatrixIndex,iItem,k,wMatrix,bVector,pSet);
v =Vi(pMatrixIndex,iItem,k,wMatrix,bVector,pSet);
return (du*v - u*dv)/(v*v);
}
实现代码在BuildDeltaW_b中
其核心就是三重循环，一是k从0-9循环，二是m从0-783循环，三是i从0到n循环，然后计算损失函数的导数。
dsm += pPi[i]*Dpi_dWmk(pMatrixIndex,i,m,k,wMatrix,bVector,pSet);
pPi[i]就是1/Pi的速查值。
test_sgd的核心就是
循环：0-Epoch
循环：0-每一个batch
计算这个batch的损失函数偏导数
然后新参数=参数-学习率*偏导数
输出新参数矩阵
如果该验证一下，就去验证一下，
如果验证的有改进，就去test一下。
判断是否可以提前结束
其他的…就不写了吧，反正test_sgd也挺简单的。
懒…
5、感谢
这部分内容对我来说是最痛苦的，特别要感谢Rachel Zhang同学。Rachel耐心、细致的回答了我无数幼稚和愚笨的基础数学问题，并鼓励我尝试完成，使得这一切成为可能。
此外，还感谢http://blog.sina.com.cn/s/blog_62339a2401015jyq.html，其图示使得我豁然开朗，明白了梯度下降算法的意义。建议初学者可以去看看。