📄 main.cpp
字号:
for(j=0;j<=nHidd;j++)
for(k=1;k<=nOut;k++)
{
pGradH[j][k]+=pNet->out_delta[k]*pNet->hidd_unit[j];
}
}
//权值调整//////////////////////////
for(j=0;j<=nIn;j++)
for(k=1;k<=nHidd;k++)
{
pNet->in_w[j][k]+=eta*pGradIn[j][k];
}
for(j=0;j<=nHidd;j++)
for(k=1;k<=nOut;k++)
{
pNet->hidd_w[j][k]+=eta*pGradH[j][k];
}
esqrsum=0.0;
//计算权值更新后的平方误差
for(p=0;p<np;p++)
//----------------------- Page 11-----------------------
{
//输入数据
pNet->in_unit[0]=0.0; //无偏移
for(i=1;i<=nIn;i++)
pNet->in_unit[i]=indata[p][i-1];
for(i=1;i<=nOut;i++)
pNet->target[i]=targetdata[p][i-1];
bpnn_nnforward(pNet);
esqr=0.0;
for(i=1;i<=nOut;i++)
{
esqr+=pNet->target[i]-pNet->out_unit[i];//计算累积平方误差
}
esqrsum+=esqr*esqr;
}
bpnn_mfree2d(pGradIn,nIn+1);
bpnn_mfree2d(pGradH,nHidd+1);
return esqrsum;
}
/**//////////////////////////////////////////////////////////////////////
//Levenberg-Marquardt的BP网络训练(全部样本)
// 参数
// pNet --已初始化的神经网络
// indata --输入的样本数组
// targetdata--输出期望
// np --样本个数
// lamda --lamda系数
// 返回
// 返回更新后的平方误差[Page]
// 更新pNet
/**//////////////////////////////////////////////////////////////////////
double bpnn_train_lm(BPNN_t *pNet,double **indata,double **targetdata,int np,double *lamda)
{
int i,j,k,p,pos;
double errout,errh;
double esqrsum,esqr;
double *pGrad; //合并的梯度向量
double **pGradIn,**pGradH; //每一层的梯度矩阵
//----------------------- Page 12-----------------------
double **pGradIn2,**pGradH2;
double *pHess; //H矩阵
double *pAlpha,*pBeta; //用于计算逆矩阵
int nIn,nHidd,nOut;
int vecsize;
nIn=pNet->in_n;
nHidd=pNet->hidd_n;
nOut=pNet->out_n;
vecsize=(nIn+1)*(nHidd)+(nHidd+1)*(nOut);
pGrad=bpnn_malloc1d(vecsize);
pHess=bpnn_malloc1d(vecsize*vecsize);
pAlpha=bpnn_malloc1d(vecsize*vecsize);
pBeta=bpnn_malloc1d(vecsize);
pGradIn=bpnn_malloc2d(nIn+1,nHidd+1);
pGradH=bpnn_malloc2d(nHidd+1,nOut+1);
pGradIn2=bpnn_malloc2d(nIn+1,nHidd+1);
pGradH2=bpnn_malloc2d(nHidd+1,nOut+1);
//每一层梯度清0
for(j=0;j<=nIn;j++)
for(k=0;k<=nHidd;k++)
{
pGradIn[j][k]=0.0;
pGradIn2[j][k]=0.0;
}
for(j=0;j<=nHidd;j++)
for(k=0;k<=nOut;k++)
{
pGradH[j][k]=0.0;
pGradH2[j][k]=0.0;
}
//计算梯度向量
for(p=0;p<np;p++)
{
//输入数据
pNet->in_unit[0]=0.0; //无偏移
for(i=1;i<=nIn;i++)
pNet->in_unit[i]=indata[p][i-1];
for(i=1;i<=nOut;i++)
pNet->target[i]=targetdata[p][i-1];
//----------------------- Page 13-----------------------
bpnn_nnforward(pNet);
bpnn_out_error(pNet->out_delta,pNet->out_dfda,pNet->target,pNet->out_unit,nOut,&errout);
bpnn_hidd_error(pNet->hidd_delta,pNet->out_delta,pNet->hidd_dfda,pNet->out_dfda,nHidd,nOut,pNet->hidd_w,pNet->hidd_unit,&errh);
pos=0;
//累加梯度
for(j=0;j<=nIn;j++)
for(k=1;k<=nHidd;k++)
{
pGradIn[j][k]+=pNet->hidd_dfda[k]*pNet->in_unit[j];
pGradIn2[j][k]+=pNet->hidd_delta[k]*pNet->in_unit[j];
pGrad[pos]=pGradIn[j][k]; //归并到梯度向量
pBeta[pos]=pGradIn2[j][k];
++pos;
}
for(j=0;j<=nHidd;j++)
for(k=1;k<=nOut;k++)
{
pGradH[j][k]+=pNet->out_dfda[k]*pNet->hidd_unit[j];
pGradH2[j][k]+=pNet->out_delta[k]*pNet->hidd_unit[j];
pGrad[pos]=pGradH[j][k]; //归并到梯度向量
pBeta[pos]=pGradH2[j][k];
++pos;
}
}//endfor
//构造H矩阵
for(j=0;j<vecsize;j++)
{
for(i=0;i<vecsize;i++)
{
pHess[i+j*vecsize]=pGrad[j]*pGrad[i];
pAlpha[i+j*vecsize]=pHess[i+j*vecsize];
}
//添加对角线
pAlpha[j+j*vecsize]=pHess[j+j*vecsize]*(1.0+*lamda);
}
//算出权值增量
gauss_jordan(pAlpha,pBeta,vecsize);
pos=0;
//----------------------- Page 14-----------------------
//权值调整//////////////////////////
for(j=0;j<=nIn;j++)
for(k=1;k<=nHidd;k++)
{
pNet->prev_in_w[j][k]=pNet->in_w[j][k];
pNet->in_w[j][k]+=pBeta[pos];
++pos;
}
for(j=0;j<=nHidd;j++)
for(k=1;k<=nOut;k++)
{
pNet->prev_hidd_w[j][k]=pNet->hidd_w[j][k];
pNet->hidd_w[j][k]+=pBeta[pos];
++pos;
}
//计算权值更新后的平方误差
esqrsum=0.0;
for(p=0;p<np;p++)
{
//输入数据
pNet->in_unit[0]=0.0; //无偏移
for(i=1;i<=nIn;i++)
pNet->in_unit[i]=indata[p][i-1];
for(i=1;i<=nOut;i++)
pNet->target[i]=targetdata[p][i-1];
bpnn_nnforward(pNet);
esqr=0.0;
for(i=1;i<=nOut;i++)
{
esqr+=pNet->target[i]-pNet->out_unit[i];//计算累积平方误差
}
esqrsum+=esqr*esqr;
}
//释放内存
bpnn_mfree2d(pGradIn,nIn+1);
bpnn_mfree2d(pGradH,nHidd+1);
bpnn_mfree2d(pGradIn2,nIn+1);
bpnn_mfree2d(pGradH2,nHidd+1);
free(pHess);
free(pBeta);
//----------------------- Page 15-----------------------
free(pAlpha);
free(pGrad);
return esqrsum;
}
/**/////////////////////////////////////////////////////////////////
// gauss-jordan消元法解线性方程组Ax=B
// 参数
// a--输入系数矩阵，输出被其逆矩阵代替(n*n)[Page]
// b--等号右边的矩阵B，输出时被对应解向量代替(n*1)
// n-- a的维数
//
/**//////////////////////////////////////////////////////////////////
void gauss_jordan(double *a,double *b,int n)
{
int i,j,k,l,ll;
int icol,irow;
double big,dum,pivinv;
int *ipiv=NULL;
double *inva=NULL;
ipiv=(int*)malloc(n*sizeof(int));
inva=(double*)malloc(n*n*sizeof(double));
for(j=0;j<n;j++)
{
ipiv[j]=0;
}
SetEye(inva,n);
for(i=0;i<n;i++)
{
big=0.0;
for(j=0;j<n;j++) //寻找主元(最大数)
{
if(ipiv[j]!=1)
for(k=0;k<n;k++)
{
if(ipiv[k]==0)
{
if(fabs(a[k+j*n])>=big)
//----------------------- Page 16-----------------------
{
big=fabs(a[k+j*n]);
irow=j;
icol=k;
}
}//endifipiv
}
}//endforj
++(ipiv[icol]);
//实施消元
if(irow!=icol)
{
SWAP(b[irow],b[icol]);
for(l=0;l<n;l++)
{
SWAP(a[l+irow*n],a[l+icol*n]);
SWAP(inva[l+irow*n],inva[l+icol*n]);
}
}
if(fabs(a[icol+icol*n])<EPSILON)
{
return; //error奇异矩阵
}
pivinv=1.0/a[icol+icol*n];
b[icol]*=pivinv;
for(l=0;l<n;l++)
{
a[l+icol*n]*=pivinv;
inva[l+icol*n]*=pivinv;
}
for(ll=0;ll<n;ll++)
{
if(ll!=icol)
{
dum=a[icol+ll*n];
b[ll]-=b[icol]*dum;
for(l=0;l<n;l++)
{
a[l+ll*n]-=a[l+icol*n]*dum;
inva[l+ll*n]-=inva[l+icol*n]*dum;
}
}
//----------------------- Page 17-----------------------
}//endfor
}//endforwhole
//更新a为逆矩阵
for(i=0;i<n*n;i++)
a[i]=inva[i];
free(ipiv);
free(inva);
}
//设置单位矩阵
void SetEye(double *mat,int n)
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i==j)
mat[j+i*n]=1.0;
else
mat[j+i*n]=0.0;
}
}
//测试程序
//testprogram
int main(int argc,char *argv[])
{
const double prec=0.001; //精度
int maxtry=5000;
int cnt=0;
int i,j,k;
double lamda=0.1;
BPNN_t *pNNet;
double **x;
double **y;
double esqr,oesqr,res;
int nIn,nHidd,nOut;
int nsample=20; //样本数
//初始化样本数据
x=bpnn_malloc2d(nsample,1);
//----------------------- Page 18-----------------------
y=bpnn_malloc2d(nsample,1);
for(i=0;i<nsample;i++)
{
x[i][0]=(double)i/(double)nsample;
//用正弦曲线验证
y[i][0]=0.5*sin(x[i][0]*PI*0.5)+0.5;
//y[i][0]=-0.4*x[i][0]+1;
//y[i][0]=0.25*(x[i][0]-2)*(x[i][0]-2);
}
//初始化神经网络
pNNet=bpnn_create(1,8,1);
nIn=pNNet->in_n;
nHidd=pNNet->hidd_n;
nOut=pNNet->out_n;
oesqr=9999; //infinite
while(cnt++<maxtry)
{
esqr=bpnn_train_lm(pNNet,x,y,nsample,&lamda);
if(esqr<prec)
break;
else
{
if(esqr<oesqr) //接受新解[Page]
{
lamda*=0.1; //缩小lamda
oesqr=esqr;
}
else
{
lamda*=10.0; //扩大lamda
//恢复旧权值
for(j=0;j<=nIn;j++)
for(k=1;k<=nHidd;k++)
{
pNNet->in_w[j][k]=pNNet->prev_in_w[j][k];
}
for(j=0;j<=nHidd;j++)
for(k=1;k<=nOut;k++)
{
pNNet->hidd_w[j][k]=pNNet->prev_hidd_w[j][k];
}
//----------------------- Page 19-----------------------
}
}
/**//*传统BP
esqr=bpnn_train_steepest(pNNet,0.1,x,y,nsample);
if(esqr<prec)
break;
*/
}
//输出测试数据
printf("输入样本:期望值 训练后的结果");
for(i=0;i<nsample;i++)
{
pNNet->in_unit[1]=x[i][0];
bpnn_nnforward(pNNet);
res=pNNet->out_unit[1];
printf("%5f : %5f %6f",x[i][0],y[i][0],res);
}
printf(" esqr:%5f ",esqr);
bpnn_destroy(pNNet);
bpnn_mfree2d(x,nsample);
bpnn_mfree2d(y,nsample);
return 0;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -