?? 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 + -