bpnet.txt

和好的算法

        //累加梯度
        for(j=0; j<=nIn; j++)
        for(k=1; k<=nHidd; k++)
        ...{
            pGradIn[j][k] += pNet->hidd_delta[k] * pNet->in_unit[j];
        }
        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++)
    ...{
        //输入数据
        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系数
//    返回
//        返回更新后的平方误差
//      更新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;    //每一层的梯度矩阵
    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];
    
        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;
        }
    }//end for

    //构造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;
    //权值调整//////////////////////////
    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);
    free(pAlpha);
    free(pGrad);

    return esqrsum;    
}



/**/////////////////////////////////////////////////////////////////
//  gauss-jordan消元法解线性方程组Ax = B
//  参数
//        a -- 输入系数矩阵，输出被其逆矩阵代替(n*n)
//      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)
                        ...{
                            big = fabs( a[k+j*n] );
                            irow = j;
                            icol = k;
                        }
                    }//end if ipiv
                }
        }//end for j
        ++(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;
                }
            }
        }//end for
    }//end for whole

    //更新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;
    }
}
 

测试程序

 

//test program
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);
    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)  //接受新解
            ...{
                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];
                }

            }
        }
        
    /**//* 传统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;
}