C 实现最小二乘,步骤以及代码

最小二乘一些复杂的公式,原理啥的我就不多说了...

直接入题: (因为符号不知道怎么打,下面的推导内容网上也找不到,所以打不出来的我全手写的)

在表达两个变量之间的关系时,有一元线性模型: Y= ax + b, 也有二次函数模型: Y = a₁X² + a₂X  + b ...等等

这里就说最简单的 Y = ax + b

最小二乘法是以实际值与目标值的差的平方和作为损失函数的,

损失函数公式: 

以上公式就是损失函数,其中的系数1/2没有实际意义,就是为了下面的推导可以把系数抹除掉.

其中的 Y = ax + b,  

x 为输入的图片特征数组,比如我上一篇博客中,Hog特征提取,详细步骤及源码 最后输出的特征大数组即可用于X的输入.

a 为输入的各个特征的权重.

a 和 x 都是多维的, b 可以看成是 a_0*1

_{其中的 T 就是目标值,这个目标值是已知的,得自己输入.}

_{上面说了 a 是各个特征对应的权重,简单点就是说一张图片里面表示的识别物体的特征的权重就会高,背景什么的权重就会低,}

_{这个权重我们也是不知道的呀,也不是和 x 一样可以有现成的值输入, 所以我们一开始就得先初始化一下这个权重 a,}

_{初始化的 a ,大家知道,保证是不正确的, 那么我们就得不断的去更新它,} 

_{更新函数:} 

_{以上公式原理我也不是很清楚,最后一个像 y 的表示学习速率,这个也是需要我们去控制的,}

_{其实我们最终的目的就是为了不断的更新权重 a ,可以使得识别物体的特征越来越明显,}

_{跟新这个 a , 我们就得求出 :
下面就是过程,详细的我就不说了,}

前面已经说过了, a,x都是多维的,所以Y也是多维的,T也是对应每个样本的目标值,得出的最后结论就是:

当然上面的更新函数也OK了,输入也都有了~

下面看看代码:

因为还没和 HOG 特征提取串联,所以输入的知识图片的imagedata, 不是特征大数组

代码很简单,就是按照上面的内容写的,最重要的就是更新权重,我用的都是opencv1.0,

样本文件夹是放在桌面,一共2套样本,一套样本目标值是-1,另一套是1

其中学习速率需要自己去设定一个合适的值,如果输出的结果不对,改下学习速率看看,

#include <iostream>#include <cstdio>#include <cv.h>#include <cxcore.h>#include <time.h>#include <highgui.h>#include <direct.h>#include <io.h>using namespace std;#define MaxLearnRate  0.00004 // 定义的学习速率, 这个数字不一定是最合适你的样本的typedef struct _SampleInfo // 每个样本对应的 data 以及目标值, 图片缩放到20*20{float iTarget;float pdImageData[400];}SampleInfo;static SampleInfo pstTrainSampleInfo[1770]; // 存放样本的数组,一共1770个样本static SampleInfo pstTestSampleInfo[3981]; // 存放样本的数组,一共1770个样本static double pdLearnRate = MaxLearnRate; // 学习速率static int iTotalTrainSample; //所有训练样本的数量static int iTotalTestSample; //所有测试样本的数量static int iTotalIteNvm = 500; // 训练的轮数static int GetTrainSumNvm(char* pSrcPath, int iFlag) // 获取所有样本个数{char pSrcSubFilePath[256];struct _finddata_t FileInfo;IplImage* pImgSrc = NULL;intptr_t lHandle;int ii, iTrainNvm=0;IplImage* pResizeImage = cvCreateImage(cvSize(20,20), 8, 1);uchar* pucImgData = (uchar*)pResizeImage->imageData;for (ii=0; ii<2; ii++)// 有 2 个文件夹装有样本{sprintf_s( pSrcSubFilePath, 256, "%s\\%d\\*.bmp", pSrcPath, ii );lHandle = _findfirst( pSrcSubFilePath, &FileInfo );if ( -1 == lHandle ){printf( "%s 路径不正确, 或者没有bmp的文件！请核实\n", pSrcSubFilePath );continue;}do{if( !( FileInfo.attrib & _A_SUBDIR ) ){sprintf_s(pSrcSubFilePath, 256, "%s\\%d\\%s", pSrcPath, ii, FileInfo.name);pImgSrc = cvLoadImage(pSrcSubFilePath, CV_LOAD_IMAGE_GRAYSCALE);cvResize(pImgSrc, pResizeImage,CV_INTER_LINEAR);if (iFlag == 1){for (int jj=0; jj<400; jj++){pstTrainSampleInfo[iTrainNvm].pdImageData[jj] = (pucImgData[jj] / 128.0f) - 1.0f;}pstTrainSampleInfo[iTrainNvm].iTarget = (ii==0) ? -1.0f : 1.0f;}else{for (int jj=0; jj<400; jj++){pstTestSampleInfo[iTrainNvm].pdImageData[jj] = (pucImgData[jj] / 128.0f) - 1.0f;}pstTestSampleInfo[iTrainNvm].iTarget = (ii==0) ? -1.0f : 1.0f;}iTrainNvm++;cvReleaseImage(&pImgSrc);}} while ( _findnext( lHandle, &FileInfo ) == 0 );}return iTrainNvm;cvReleaseImage(&pResizeImage);}int GradientDescend() {CvMat* cmTrainData = cvCreateMat(iTotalTrainSample, 401, CV_32FC1); // 401 是因为 还有一个 a0*1, 即 Y = ax + b 中的 bCvMat* cmTrainResult = cvCreateMat(iTotalTrainSample, 1, CV_32FC1);CvMat* cmTrainTarget = cvCreateMat(iTotalTrainSample, 1, CV_32FC1);CvMat* cmSubVal = cvCreateMat(iTotalTrainSample, 1, CV_32FC1);CvMat* cmSubValT = cvCreateMat(1, iTotalTrainSample, CV_32FC1);CvMat* cmWeight = cvCreateMat(1, 401, CV_32FC1);CvMat* cmWeightT = cvCreateMat(401, 1, CV_32FC1);CvMat* cmDataMulSub = cvCreateMat(1, 401, CV_32FC1);// 初始化权重for(int ii = 0; ii < 401; ii++){cmWeight->data.fl[ii]= ( rand() % 1000 ) / 1000.0 - 0.5;}// 将样本中的数据转化到矩阵中for (int mm=0; mm<iTotalTrainSample; mm++){memcpy_s(cmTrainData->data.fl + mm * 401, sizeof(float)*400, pstTrainSampleInfo[mm].pdImageData, sizeof(float)*400);cmTrainTarget->data.fl[mm] = pstTrainSampleInfo[mm].iTarget;cmTrainData->data.fl[mm*401 + 400] = 1.0f; // 将所有样本中的 b 赋值成 1}int iPosError = 0, iNegError = 0, iLastErrorSum=0; cvTranspose(cmWeight, cmWeightT);for (int ii = 0; ii < iTotalIteNvm; ii++) {iPosError = iNegError = 0;//// 下面的步骤就是更新权重,cvMatMul(cmTrainData, cmWeightT, cmTrainResult); // Y= ax + b , cmData中 已经包含了 b(a0*1),cvSub(cmTrainResult, cmTrainTarget, cmSubVal, NULL); // Y - T , 实际值 - 目标值cvTranspose(cmSubVal, cmSubValT);cvMatMul(cmSubValT, cmTrainData, cmDataMulSub); // (Y- T) * xcvConvertScale(cmDataMulSub, cmDataMulSub, pdLearnRate, 0);// (Y- T) * x * 学习速率cvSub(cmWeight, cmDataMulSub, cmWeight, NULL); // 更新权重 cvTranspose(cmWeight, cmWeightT);cvMatMul(cmTrainData, cmWeightT, cmTrainResult); // 权重更新完成后, 重新 Y= ax + b// 求出的实际值 和 目标值进行比较, 得出训练集的错误样本的个数for (int mm=0; mm<iTotalTrainSample; mm++){if ( cmTrainResult->data.fl[mm] < 0 && cmTrainTarget->data.fl[mm] == 1 ){iPosError++;}if ( cmTrainResult->data.fl[mm]  > 0 && cmTrainTarget->data.fl[mm] == -1 ){iNegError++;}}printf("train ------  %d  -> %d  %d  %d  -> %f \n", iTotalTrainSample, iPosError, iNegError, (iPosError+iNegError) , ((float)(iPosError+iNegError))/(float)iTotalTrainSample);//这一轮训练的误差比上一轮大，此时应该降低学习速率。if ( ii!=0 && (iPosError+iNegError)>=iLastErrorSum){pdLearnRate = (pdLearnRate<MaxLearnRate/100000) ? MaxLearnRate: pdLearnRate / 2;}iLastErrorSum = iPosError + iNegError;}printf("============================================================\n" );// 下面是根据上面求出的 cmWeightT(权重) 得出测试集的错误个数CvMat* cmTestData = cvCreateMat(iTotalTestSample, 401, CV_32FC1); // 401 是因为 还有一个 a0*1, 即 Y = ax + b 中的 bCvMat* cmTestResult = cvCreateMat(iTotalTestSample, 1, CV_32FC1);CvMat* cmTestTarget = cvCreateMat(iTotalTestSample, 1, CV_32FC1);int iPosError_t = 0, iNegError_t = 0;//, iLastErrorSum_t =0; for (int mm=0; mm<iTotalTestSample; mm++){memcpy_s(cmTestData->data.fl + mm * 401, sizeof(float)*400, pstTestSampleInfo[mm].pdImageData, sizeof(float)*400);cmTestTarget->data.fl[mm] = pstTestSampleInfo[mm].iTarget;cmTestData->data.fl[mm*401 + 400] = 1.0f; // 将所有样本中的 b 赋值成 1}iPosError_t = iNegError_t = 0;cvMatMul(cmTestData, cmWeightT, cmTestResult); // 权重更新完成后, 重新 Y= ax + b// 求出的实际值 和 目标值进行比较, 得出错误样本的个数for (int mm=0; mm<iTotalTestSample; mm++){if ( cmTestResult->data.fl[mm] < 0 && cmTestTarget->data.fl[mm] == 1 ){iPosError_t++;}if ( cmTestResult->data.fl[mm]  > 0 && cmTestTarget->data.fl[mm] == -1 ){iNegError_t++;}} printf("Test ------  %d  ->  %d  %d  %d -> %f \n",iTotalTestSample, iPosError_t, iNegError_t, (iPosError_t+iNegError_t), ((float)(iPosError_t+iNegError_t))/(float)iTotalTestSample );return 0;}int main(){// 样本放在桌面iTotalTrainSample = GetTrainSumNvm("C:\\Users\\Administrator\\Desktop\\Train", 1);iTotalTestSample = GetTrainSumNvm("C:\\Users\\Administrator\\Desktop\\Test", 0);srand((int)time(0));GradientDescend();system("pause");return 0;}

下面是结果:

因为输入的直接是 imagedata 所以错误率有点高 将近12%了,

如果输入为 上篇博客中 的 Hog特征数组,那么错误率就会小很多~

没事干吧Hog特征和最小二乘整合了一个分类器,代码就是上一篇博客和这篇博客合起来的,我就不贴了,

贴下效果图,错误率是小了很多,但是还不够小~所以继续学习~~

这次是4个文件夹的样本进行训练和测试,每个样本的权重是4维的,属于第一类就是[1,-1,-1,-1],

这样计算出来有四个实际值,取其中的最大值和目标值进行比较来确定是不是错误样本,

训练迭代了2000次的效果图:

在以上基础上,我加了个 双曲正切激励激励函数 tanh(), 迭代到1000次的时候就有上图的结果,

下面是迭代2000的结果:

我觉得还有戏,直接迭代到1W次,效果不错呦~Hog特征还是比较猛的~

迭代到1W次的结果: