nnssif.c
来自「基于MATLAB的神经网络非线性系统辨识软件包.」· C语言 代码 · 共 818 行 · 第 1/3 页
C
818 行
/*
>>>>>>>>>>>>>>>>>>>>>>>>>>> COMPUTE THE PSI MATRIX <<<<<<<<<<<<<<<<<<<<<<<<<<
(The derivative of each network output (y2) with respect to each weight)
*/
/* Some intermidiate computations */
for(j=0;j<hhids;j++)
{
jj = (int)cvget(H_hidden,j);
for(k=0;k<N;k++)
put_val(tmp0,j,k,1-get_val(y1,jj,k)*get_val(y1,jj,k));
}
/* ========== Elements corresponding to the linear output units ===========*/
for(i=0; i<louts; i++)
{
ii = (int)cvget(L_output,i);
/*** The part of PSIx corresponding to hidden-to-output layer weights ***/
index1 = ii * (hidden+1);
psi1(PSIx, index1, index2, ii, y1);
/************************************************************************/
/**** The part of PSIx corresponding to input-to-hidden layer weights ****/
for(j=0; j<lhids; j++)
{
jj = (int)cvget(L_hidden,j);
psi2(PSIx, (int)cvget(index,jj), index2, ii, get_val(W2,ii,jj), PHI);
}
for(j=0; j<hhids;j++)
{
jj = (int)cvget(H_hidden,j);
psi3(tmp3, tmp0, j, get_val(W2,ii,jj));
psi4(PSIx, (int)cvget(index,jj), index2, ii, tmp3, PHI);
}
/************************************************************************/
}
/* =========== Elements corresponding to the tanh output units ===========*/
for(i=0; i<houts; i++)
{
ii = (int)cvget(H_output,i);
index1 = ii * (hidden + 1);
for(k=0; k<N; k++)
put_val(tmp2,0,k,1-get_val(y2,ii,k)*get_val(y2,ii,k));
/* -- The part of PSIx corresponding to hidden-to-output layer weights --*/
psi4(PSIx, index1, index2, ii, tmp2, y1);
/* ---------------------------------------------------------------------*/
/* -- The part of PSIx corresponding to input-to-hidden layer weights ---*/
for(j=0; j<lhids; j++)
{
jj = (int)cvget(L_hidden,j);
smul(tmp3, tmp2, get_val(W2,ii,jj));
psi4(PSIx, (int)cvget(index,jj), index2, ii, tmp3, PHI);
}
for(j=0; j<hhids; j++)
{
jj = (int)cvget(H_hidden,j);
psi3(tmp3, tmp0, j, get_val(W2,ii,jj));
psi5(PSIx, (int)cvget(index,jj), index2, ii, tmp3, tmp2, PHI);
}
/* ---------------------------------------------------------------------*/
}
/*
>>>>>>>>>>>>>>>>>>>> Linearize network <<<<<<<<<<<<<<<<<<<<<
*/
for(t=0;t<N;t++){
/*-- Derivative of states wrt. hidden outputs --*/
for(j=0;j<louts;j++){
i=(int)cvget(L_output,j);
for(k=0;k<hidden;k++) put_val(dxdy1,i,k,get_val(W2,i,k));
}
for(j=0;j<houts;j++){
i=(int)cvget(H_output,j);
for(k=0;k<hidden;k++) put_val(dxdy1,i,k,get_val(W2,i,k)*(1-\
get_val(y2,i,t)*get_val(y2,i,t)));
}
/*-- Partial deriv. of output from each hidden unit wrt. net inputs --*/
for(j=0;j<lhids;j++){
i=(int)cvget(L_hidden,j);
for(k=0;k<nx;k++) put_val(dy1dx,i,k,get_val(W1,i,k));
for(k=nxu;k<inputs;k++) put_val(dy1de,i,k-nxu,get_val(W1,i,k));
}
for(j=0;j<hhids;j++){
i=(int)cvget(H_hidden,j);
for(k=0;k<nx;k++) put_val(dy1dx,i,k,\
get_val(W1,i,k)*(1-get_val(y1,i,t)*get_val(y1,i,t)));
for(k=nxu;k<inputs;k++) put_val(dy1de,i,\
k-nxu,get_val(W1,i,k)*(1-get_val(y1,i,t)*get_val(y1,i,t)));
}
/*--Partial derivative of states w.r.t. past states and residuals --*/
mmul(Ahat,dxdy1,dy1dx);
for(k=0;k<nx-ny;k++) put_val(Ahat,(int)cvget(nrowidx,k),\
(int)cvget(nrowidx,k)+1,1.0);
mmul(Khat,dxdy1,dy1de);
mmul(AKC,Khat,C);
msub(AKC,Ahat,AKC);
/*
>>>>>>>>>>>>>>>>>>> Filter partial derivatives <<<<<<<<<<<<<<<<<<<<
*/
if(t>=1){
/* PSIx = PSIx + PSIx*AKC' */
index5 = t*nx;
index6 = (t-1)*nx;
for(i=0;i<reduced;i++){
ii =(int)cvget(theta_index,i);
for(j=0;j<nx;j++){
for(k=0;k<nx;k++){
PSIx->mat[ii][index5+j]+=get_val(PSIx,ii,index6+k)*\
get_val(AKC,j,k);
}
}
}
}
/*PSI=PSIx*C';*/
index5 = t*ny;
index6 = t*nx;
for(i=0;i<reduced;i++){
ii =(int)cvget(theta_index,i);
for(j=0;j<ny;j++){
for(sum=0,k=0;k<nx;k++){
sum+=get_val(PSIx,ii,index6+k)*get_val(C,j,k);
}
put_val(PSI,ii,index5+j,sum);
}
}
}
minit(PSIx);
dw = 0;
/*
>>>>>>>>>>>> Gradient (G = PSI_red*E_vector - D*theta_red) <<<<<<<<<<<<<
*/
for(i=0; i<reduced; i++){
ii = (int)cvget(theta_index,i);
for(sum=0.0,k=skipstart; k<Nny; k++)
sum+=get_val(PSI,ii,k)*rvget(Evec,k);
cvput(G,i,sum - cvget(D,i)*cvget(theta_red,i));
}
/*
>>>>>>>>>> Mean square error part of Hessian (PSI_red*PSI_red') <<<<<<<<<<
*/
for(i=0; i<reduced; i++){
ii = (int)cvget(theta_index,i);
for(j=i; j<reduced; j++){
jj = (int)cvget(theta_index,j);
for(sum=0.0,k=skipstart; k<Nny; k++)
sum += get_val(PSI,ii,k)*get_val(PSI,jj,k);
put_val(H,i,j,sum);
put_val(H,j,i,sum);
}
}
for(i=0;i<reduced;i++) /* Add diagonal matrix */
put_val(H,i,i,get_val(H,i,i)+cvget(D,i));
}
/*
>>>>>>>>>>>>>>>>>>>>>>>>>>> COMPUTE h_k <<<<<<<<<<<<<<<<<<<<<<<<<<<
*/
/* -- Hessian (H = R + lambda*I + D) --*/
tmp1 = lambda - lambda_old;
for(i=0;i<reduced;i++) /* Add diagonal matrix */
put_val(H,i,i,get_val(H,i,i)+tmp1);
/* -- Search direction -- */
choldc(H, Htmp);
cholsl(Htmp,h,G);
/* -- Compute 'apriori' iterate -- */
madd(theta_red_new,theta_red,h); /* Update parameter vector */
mcopyi(theta,theta_index,index0,theta_red_new,index7,index0);
/* -- Put the parameters back into the weight matrices -- */
v2mreshape(W1_new,theta,parameters2);
v2mreshape(W2_new,theta,0);
/*
>>>>>>>>>>>>> Compute network output y2(theta+h) <<<<<<<<<<<<<<
*/
for(t=0;t<N;t++){
mvmul(h1,W1_new,PHI,t);
vtanh(y1,H_hidden,t,h1,H_hidden,0);
vcopyi(y1,L_hidden,t,h1,L_hidden,0);
mvmul(h2,W2_new,y1,t);
vtanh(y2,H_output,t,h2,H_output,0);
vcopyi(y2,L_output,t,h2,L_output,0);
for(k=0;k<(nx-ny);k++){
i=(int)cvget(nrowidx,k);
y2->mat[i][t]+=get_val(PHI,i+1,t);
}
mvmul(Yhat,C,y2,t); /* Output prediction */
for(i=0;i<ny;i++){
cvput(E,i,get_val(Y2,i,t)-cvget(Yhat,i));/* Prediction error */
rvput(Evec_new,t*ny+i,cvget(E,i)); /* Store E in Evec_new */
}
if(t<N-1){
for(i=0;i<nx;i++)
put_val(PHI,i,t+1,get_val(y2,i,t));
for(i=0;i<ny;i++)
put_val(PHI,nx+nu+i,t+1,cvget(E,i));
}
}
for(SSE_new=0,t=skipstart;t<Nny;t++)
SSE_new+=rvget(Evec_new,t)*rvget(Evec_new,t); /* Sum of squared errors */
for(tmp1=0,i=0;i<reduced;i++) tmp1+=cvget(theta_red_new,i)*cvget(theta_red_new,i)*cvget(D,i);
NSSE_new = (SSE_new+tmp1)/(2*N2); /* Value of cost function*/
/*
>>>>>>>>>>>>>>>>>>>>>>>>>>> UPDATE lambda <<<<<<<<<<<<<<<<<<<<<<<<<<
*/
lambda_old = lambda;
for(tmp1=0,i=0;i<reduced;i++) tmp1+=cvget(h,i)*cvget(h,i)*(cvget(D,i)+lambda);
L = sprod3(h,G) + tmp1;
/* Decrease lambda if SSE has fallen 'sufficiently' */
if(2*N2*(NSSE - NSSE_new) > (0.75*L)) lambda = lambda/2;
/* Increase lambda if SSE has grown 'sufficiently' */
else if(2*N2*(NSSE-NSSE_new) <= (0.25*L)) lambda = 2*lambda;
/*
>>>>>>>>>>>>>>>>>>> UPDATES FOR NEXT ITERATION <<<<<<<<<<<<<<<<<<<<
*/
/* Update only if criterion has decreased */
if(NSSE_new<NSSE)
{
critdif = NSSE-NSSE_new; /* Criterion difference */
for(i=0,gradmax=0.0,ptm1=G->mat[0];i<reduced;i++){ /* Maximum gradient */
sum = fabs(*(ptm1++));
if(gradmax<sum)
gradmax = sum;
}
gradmax/=N2;
ptm1=theta_red_new->mat[0];
ptm2=theta_red->mat[0];
for(i=0,paramdif=0.0;i<reduced;i++){ /* Maximum gradient */
sum = fabs(*(ptm1++) - *(ptm2++));
if(paramdif<sum)
paramdif = sum;
}
lambda_old = 0.0;
tmp = W1; W1=W1_new; W1_new=tmp;
tmp = W2; W2=W2_new; W2_new=tmp;
tmp = theta_red; theta_red=theta_red_new; theta_red_new = tmp;
tmp = Evec; Evec = Evec_new; Evec_new = tmp;
dw = 1;
NSSE = NSSE_new;
cvput(NSSEvec,iteration-1,NSSE);
switch(trparms->infolevel){ /* Print on-line inform */
case 1:
printf("# %i W=%4.3e critdif=%3.2e maxgrad=%3.2e paramdif=%3.2e\n",
iteration,NSSE,critdif,gradmax,paramdif);
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?