径向基函数（RBF）神经网络

径向基函数（RBF）神经网络

RBF网络能够逼近任意的非线性函数，可以处理系统内的难以解析的规律性，具有良好的泛化能力，并有很快的学习收敛速度，已成功应用于非线性函数逼近、时间序列分析、数据分类、模式识别、信息处理、图像处理、系统建模、控制和故障诊断等。

简单说明一下为什么RBF网络学习收敛得比较快。当网络的一个或多个可调参数（权值或阈值）对任何一个输出都有影响时，这样的网络称为全局逼近网络。由于对于每次输入，网络上的每一个权值都要调整，从而导致全局逼近网络的学习速度很慢。BP网络就是一个典型的例子。

如果对于输入空间的某个局部区域只有少数几个连接权值影响输出，则该网络称为局部逼近网络。常见的局部逼近网络有RBF网络、小脑模型（CMAC）网络、B样条网络等。

径向基函数解决插值问题

完全内插法要求插值函数经过每个样本点，即。样本点总共有P个。

RBF的方法是要选择P个基函数，每个基函数对应一个训练数据，各基函数形式为，由于距离是径向同性的，因此称为径向基函数。||X-X^p||表示差向量的模，或者叫2范数。

基于为径向基函数的插值函数为：

输入X是个m维的向量，样本容量为P，P>m。可以看到输入数据点X^p是径向基函数φ_p的中心。

隐藏层的作用是把向量从低维m映射到高维P，低维线性不可分的情况到高维就线性可分了。

将插值条件代入：

写成向量的形式为，显然Φ是个规模这P对称矩阵，且与X的维度无关，当Φ可逆时，有。

对于一大类函数，当输入的X各不相同时，Φ就是可逆的。下面的几个函数就属于这“一大类”函数：

1）Gauss（高斯）函数

2）Reflected Sigmoidal（反常S型）函数

3）Inverse multiquadrics（拟多二次）函数

σ称为径向基函数的扩展常数，它反应了函数图像的宽度，σ越小，宽度越窄，函数越具有选择性。

完全内插存在一些问题：

1）插值曲面必须经过所有样本点，当样本中包含噪声时，神经网络将拟合出一个错误的曲面，从而使泛化能力下降。

由于输入样本中包含噪声，所以我们可以设计隐藏层大小为K，K<P，从样本中选取K个（假设不包含噪声）作为Φ函数的中心。

2）基函数个数等于训练样本数目，当训练样本数远远大于物理过程中固有的自由度时，问题就称为超定的，插值矩阵求逆时可能导致不稳定。

拟合函数F的重建问题满足以下3个条件时，称问题为适定的：

1. 解的存在性
2. 解的唯一性
3. 解的连续性

不适定问题大量存在，为解决这个问题，就引入了正则化理论。

正则化理论

正则化的基本思想是通过加入一个含有解的先验知识的约束来控制映射函数的光滑性，这样相似的输入就对应着相似的输出。

寻找逼近函数F(x)通过最小化下面的目标函数来实现：

加式的第一项好理解，这是均方误差，寻找最优的逼近函数，自然要使均方误差最小。第二项是用来控制逼近函数光滑程度的，称为正则化项，λ是正则化参数，D是一个线性微分算子，代表了对F(x)的先验知识。曲率过大（光滑度过低）的F(x)通常具有较大的||DF||值，因此将受到较大的惩罚。

直接给出(1)式的解：

权向量********************************(2)

G(X,X^p)称为Green函数，G称为Green矩阵。Green函数与算子D的形式有关，当D具有旋转不变性和平移不变性时，。这类Green函数的一个重要例子是多元Gauss函数：

。

正则化RBF网络

输入样本有P个时，隐藏层神经元数目为P，且第p个神经元采用的变换函数为G(X,X^p)，它们相同的扩展常数σ。输出层神经元直接把净输入作为输出。输入层到隐藏层的权值全设为1,隐藏层到输出层的权值是需要训练得到的：逐一输入所有的样本，计算隐藏层上所有的Green函数，根据(2)式计算权值。

广义RBF网络

Cover定理指出：将复杂的模式分类问题非线性地映射到高维空间将比投影到低维空间更可能线性可分。

广义RBF网络：从输入层到隐藏层相当于是把低维空间的数据映射到高维空间，输入层细胞个数为样本的维度，所以隐藏层细胞个数一定要比输入层细胞个数多。从隐藏层到输出层是对高维空间的数据进行线性分类的过程，可以采用单层感知器常用的那些学习规则，参见神经网络基础和感知器。

注意广义RBF网络只要求隐藏层神经元个数大于输入层神经元个数，并没有要求等于输入样本个数，实际上它比样本数目要少得多。因为在标准RBF网络中，当样本数目很大时，就需要很多基函数，权值矩阵就会很大，计算复杂且容易产生病态问题。另外广RBF网与传统RBF网相比，还有以下不同：

1. 径向基函数的中心不再限制在输入数据点上，而由训练算法确定。
2. 各径向基函数的扩展常数不再统一，而由训练算法确定。
3. 输出函数的线性变换中包含阈值参数，用于补偿基函数在样本集上的平均值与目标值之间的差别。

因此广义RBF网络的设计包括：

结构设计--隐藏层含有几个节点合适

参数设计--各基函数的数据中心及扩展常数、输出节点的权值。

下面给出计算数据中心的两种方法：

1. 数据中心从样本中选取。样本密集的地方多采集一些。各基函数采用统一的偏扩展常数：
   
   
   d_max是所选数据中心之间的最大距离，M是数据中心的个数。扩展常数这么计算是为了避免径向基函数太尖或太平。
2. 自组织选择法，比如对样本进行聚类、梯度训练法、资源分配网络等。各聚类中心确定以后，根据各中心之间的距离确定对应径向基函数的扩展常数。

   
   
   λ是重叠系数。

接下来求权值W时就不能再用了，因为对于广义RBF网络，其行数大于列数，此时可以求Φ伪逆。

 数据中心的监督学习算法

最一般的情况，RBF函数中心、扩展常数、输出权值都应该采用监督学习算法进行训练，经历一个误差修正学习的过程，与BP网络的学习原理一样。同样采用梯度下降法，定义目标函数为

e_i为输入第i个样本时的误差信号。

上式的输出函数中忽略了阈值。

为使目标函数最小化，各参数的修正量应与其负梯度成正比，即

具体计算式为

上述目标函数是所有训练样本引起的误差总和，导出的参数修正公式是一种批处理式调整，即所有样本输入一轮后调整一次。目标函数也可以为瞬时值形式，即当前输入引起的误差

此时参数的修正值为：

---

下面我们就分别用本文最后提到的聚类的方法和数据中心的监督学习方法做一道练习题。

考虑Hermit多项式的逼近问题

训练样本这样产生：样本数P=100，x_i且服从[-4,4]上的均匀分布，样本输出为F(x_i)+e_i，e_i为添加的噪声，服从均值为0，标准差为0.1的正态分布。

（1）用聚类方法求数据中心和扩展常数，输出权值和阈值用伪逆法求解。隐藏节点数M=10，隐藏节点重叠系数λ=1，初始聚类中心取前10个训练样本。

#include<iostream>

#include<algorithm>

#include<limits>

#include<cassert>

#include<cmath>

#include<ctime>

#include<cstdlib>

#include<vector>

#include<iomanip>

#include"matrix.h"

 

using namespace std;

 

const int P=100;        //输入样本的数量

vector<double> X(P);  //输入样本

Matrix<double> Y(P,1);        //输入样本对应的期望输出

const int M=10;         //隐藏层节点数目

vector<double> center(M);       //M个Green函数的数据中心

vector<double> delta(M);        //M个Green函数的扩展常数

Matrix<double> Green(P,M);         //Green矩阵

Matrix<double> Weight(M,1);       //权值矩阵

 

/*Hermit多项式函数*/

inline double Hermit(double x){

    return 1.1*(1-x+2*x*x)*exp(-1*x*x/2);

}

 

/*产生指定区间上均匀分布的随机数*/

inline double uniform(double floor,double ceil){

    return floor+1.0*rand()/RAND_MAX*(ceil-floor);

}

 

/*产生区间[floor,ceil]上服从正态分布N[mu,sigma]的随机数*/

inline double RandomNorm(double mu,double sigma,double floor,double ceil){

    double x,prob,y;

    do{

        x=uniform(floor,ceil);

        prob=1/sqrt(2*M_PI*sigma)*exp(-1*(x-mu)*(x-mu)/(2*sigma*sigma));

        y=1.0*rand()/RAND_MAX;

    }while(y>prob);

    return x;

}

 

/*产生输入样本*/

void generateSample(){

    for(int i=0;i<P;++i){

        double in=uniform(-4,4);

        X[i]=in;

        Y.put(i,0,Hermit(in)+RandomNorm(0,0.1,-0.3,0.3));

    }

}

 

/*寻找样本离哪个中心最近*/

int nearest(const vector<double>& center,double sample){

    int rect=-1;

    double dist=numeric_limits<double>::max();

    for(int i=0;i<center.size();++i){

        if(fabs(sample-center[i])<dist){

            dist=fabs(sample-center[i]);

            rect=i;

        }

    }

    return rect;

}

 

/*计算簇的质心*/

double calCenter(const vector<double> &g){

    int len=g.size();

    double sum=0.0;

    for(int i=0;i<len;++i)

        sum+=g[i];

    return sum/len;

}

 

/*KMeans聚类法产生数据中心*/

void KMeans(){

    assert(P%M==0);

    vector<vector<double> > group(M);          //记录各个聚类中包含哪些样本

    double gap=0.001;       //聚类中心的改变量小于为个值时，迭代终止

    for(int i=0;i<M;++i){   //从P个输入样本中随机选P个作为初始聚类中心

        center[i]=X[10*i+3];     //输入是均匀分布的，所以我们均匀地选取

    }

    while(1){

        for(int i=0;i<M;++i)

            group[i].clear();   //先清空聚类信息

        for(int i=0;i<P;++i){       //把所有输入样本归到对应的簇

            int c=nearest(center,X[i]);

            group[c].push_back(X[i]);

        }

        vector<double> new_center(M);       //存储新的簇心

        for(int i=0;i<M;++i){

            vector<double> g=group[i];

            new_center[i]=calCenter(g);

        }

        bool flag=false;

        for(int i=0;i<M;++i){       //检查前后两次质心的改变量是否都小于gap

            if(fabs(new_center[i]-center[i])>gap){

                flag=true;

                break;

            }

        }

        center=new_center;

        if(!flag)

            break;

    }

}

 

/*生成Green矩阵*/

void calGreen(){

    for(int i=0;i<P;++i){

        for(int j=0;j<M;++j){

            Green.put(i,j,exp(-1.0*(X[i]-center[j])*(X[i]-center[j])/(2*delta[j]*delta[j])));

        }

    }

}

 

/*求一个矩阵的伪逆*/

Matrix<double> getGereralizedInverse(const Matrix<double> &matrix){

    return (matrix.getTranspose()*matrix).getInverse()*(matrix.getTranspose());

}

 

/*利用已训练好的神经网络，由输入得到输出*/

double getOutput(double x){

    double y=0.0;

    for(int i=0;i<M;++i)

        y+=Weight.get(i,0)*exp(-1.0*(x-center[i])*(x-center[i])/(2*delta[i]*delta[i]));

    return y;

}

 

int main(int argc,char *argv[]){<br>　　 srand(time(0));

    generateSample();       //产生输入和对应的期望输出样本

    KMeans();           //对输入进行聚类，产生聚类中心

    sort(center.begin(),center.end());      //对聚类中心（一维数据）进行排序

     

    //根据聚类中心间的距离，计算各扩展常数

    delta[0]=center[1]-center[0];      

    delta[M-1]=center[M-1]-center[M-2];

    for(int i=1;i<M-1;++i){

        double d1=center[i]-center[i-1];

        double d2=center[i+1]-center[i];

        delta[i]=d1<d2?d1:d2;

    }

     

    calGreen();     //计算Green矩阵

    Weight=getGereralizedInverse(Green)*Y;      //计算权值矩阵

     

    //根据已训练好的神经网络作几组测试

    for(int x=-4;x<5;++x){

        cout<<x<<" ";

        cout<<setprecision(8)<<setiosflags(ios::left)<<setw(15);

        cout<<getOutput(x)<<Hermit(x)<<endl;      //先输出我们预测的值，再输出真实值

    }

    return 0;

}

并且我将其中其中的C++代码改写成了M文件

%%  
% 根据以下链接中的思想，把C++代码改写成M文件  
% http://www.cnblogs.com/zhangchaoyang/articles/2591663.html  
clear;clc;  
P=101;%训练样本共P个  
X=[]; %训练输入  
Y=[]; %训练输出  
M=10; %数据中心的个数（或说隐藏层的个数）  
  
  
centers=[];%存储数据中心（或说核函数的个数）  
deltas=[]; %存储核函数的标准差  
weights=[];%存放网络的权值（或说每个核的权值）  
set = {};  %存放不同簇所包含的所有样例  
gap=0.1;   %这是用k_means法进行聚类的时候的停止规则  
%**************************************************************************  
%构造训练样本X,Y  
X=[-4:0.08:4];  
for i=1:P  
    Y(i)=1.1*(1-X(i)+2*X(i)^2)*exp(-X(i)^2/2);  
end  
Y=Y+0.1*randn(1,P);  
%%  
%**************************************************************************  
%对输入进行聚类，(获得核函数的中心)  
for i=1:M   
    %因为我们的X是均匀分布,所以初始化也为均匀的  
    centers(i)= X( i*floor( P/10 ) );  
end  
done=0;  
while(~done)  
    for i=1:M  
       set{i}=[];  
    end  
    %计算P中每个点所属的簇  
    for i=1:P  
        distance=100;%设置一个比较大的值  
        for j=1:M  
            curr=abs(X(i)-centers(j));  
            if curr<distance  
                sets=j;  
                distance=curr;  
            end  
        end  
        set{sets}=[set{sets},X(i)];%把新分类的样例添到相应的簇中  
    end  
    %重新计算每个簇的质心  
    for i=1:M  
        new_centers(i)=sum(set{i})/length(set{i});  
    end  
    %根据各簇中心的更新情况决定是否已完成循环  
    done=0;  
%     abs(centers-new_centers)  
    for i=1:M  
        if abs(centers(i)-new_centers(i))>gap  
            done=0;  
            break;  
        else  
            done=1;  
        end  
    end  
    centers=new_centers;  
end  
%计算出每个高斯核函数的标准差（重叠系数=1）  
for i=1:M  
    curr=abs( centers-centers(i) );  
    [curr_2,b]=min(curr);  
    curr(b)=100;  
    curr_2=min(curr);  
    deltas(i)=1*curr_2;  
end  
%**************************************************************************  
%根据d=sum(K*W)  
%首先构造K为P×M的  
for i=1:P  
    for j=1:M  
        curr=abs(X(i)-centers(j));  
        K(i,j)=exp( -curr^2/(2*deltas(j)^2) );  
    end  
end  
%计算权值矩阵  
weights=inv(K'*K)*K'*Y';  
%**************************************************************************  
%测试计算出函数的情况  
x_test=[-4:0.1:4];  
for i=1:length(x_test)  
    sum=0;  
    for j=1:M  
        curr=weights(j)*exp(-abs(x_test(i)-centers(j))^2/(2*deltas(j)^2));  
        sum=sum+curr;  
    end  
    y_test(i)=sum;  
end  
figure(1)  
scatter(X,Y,'k+');  
hold on;  
plot(x_test,y_test,'r.-')  

（2）用梯度下降法训练RBF网络，设η=0.001，M=10，初始权值为[-0.1,0.1]内的随机数，初始数据中心为[-4,4]内的随机数，初始扩展常数取[0.1,0.3]内的随机数，目标误差为0.9，最大训练次数为5000。

#include<iostream>

#include<cassert>

#include<cmath>

#include<ctime>

#include<cstdlib>

#include<vector>

#include<iomanip>

 

using namespace std;

 

const int P=100;        //输入样本的数量

vector<double> X(P);  //输入样本

vector<double> Y(P);      //输入样本对应的期望输出

const int M=10;         //隐藏层节点数目

vector<double> center(M);       //M个Green函数的数据中心

vector<double> delta(M);        //M个Green函数的扩展常数

double Green[P][M];         //Green矩阵

vector<double> Weight(M);       //权值矩阵

const double eta=0.001;     //学习率

const double ERR=0.9;       //目标误差

const int ITERATION_CEIL=5000;      //最大训练次数

vector<double> error(P);  //单个样本引起的误差

 

/*Hermit多项式函数*/

inline double Hermit(double x){

    return 1.1*(1-x+2*x*x)*exp(-1*x*x/2);

}

 

/*产生指定区间上均匀分布的随机数*/

inline double uniform(double floor,double ceil){

    return floor+1.0*rand()/RAND_MAX*(ceil-floor);

}

 

/*产生区间[floor,ceil]上服从正态分布N[mu,sigma]的随机数*/

inline double RandomNorm(double mu,double sigma,double floor,double ceil){

    double x,prob,y;

    do{

        x=uniform(floor,ceil);

        prob=1/sqrt(2*M_PI*sigma)*exp(-1*(x-mu)*(x-mu)/(2*sigma*sigma));

        y=1.0*rand()/RAND_MAX;

    }while(y>prob);

    return x;

}

 

/*产生输入样本*/

void generateSample(){

    for(int i=0;i<P;++i){

        double in=uniform(-4,4);

        X[i]=in;

        Y[i]=Hermit(in)+RandomNorm(0,0.1,-0.3,0.3);

    }

}

 

/*给向量赋予[floor,ceil]上的随机值*/

void initVector(vector<double> &vec,double floor,double ceil){

    for(int i=0;i<vec.size();++i)

        vec[i]=uniform(floor,ceil);

}

 

/*根据网络，由输入得到输出*/

double getOutput(double x){

    double y=0.0;

    for(int i=0;i<M;++i)

        y+=Weight[i]*exp(-1.0*(x-center[i])*(x-center[i])/(2*delta[i]*delta[i]));

    return y;

}

 

/*计算单个样本引起的误差*/

double calSingleError(int index){

    double output=getOutput(X[index]);

    return Y[index]-output;

}

 

/*计算所有训练样本引起的总误差*/

double calTotalError(){

    double rect=0.0;

    for(int i=0;i<P;++i){

        error[i]=calSingleError(i);

        rect+=error[i]*error[i];

    }

    return rect/2;

}

 

/*更新网络参数*/

void updateParam(){

    for(int j=0;j<M;++j){

        double delta_center=0.0,delta_delta=0.0,delta_weight=0.0;

        double sum1=0.0,sum2=0.0,sum3=0.0;

        for(int i=0;i<P;++i){

            sum1+=error[i]*exp(-1.0*(X[i]-center[j])*(X[i]-center[j])/(2*delta[j]*delta[j]))*(X[i]-center[j]);

            sum2+=error[i]*exp(-1.0*(X[i]-center[j])*(X[i]-center[j])/(2*delta[j]*delta[j]))*(X[i]-center[j])*(X[i]-center[j]);

            sum3+=error[i]*exp(-1.0*(X[i]-center[j])*(X[i]-center[j])/(2*delta[j]*delta[j]));

        }

        delta_center=eta*Weight[j]/(delta[j]*delta[j])*sum1;

        delta_delta=eta*Weight[j]/pow(delta[j],3)*sum2;

        delta_weight=eta*sum3;

        center[j]+=delta_center;

        delta[j]+=delta_delta;

        Weight[j]+=delta_weight;

    }

}

 

int main(int argc,char *argv[]){

    srand(time(0));

    /*初始化网络参数*/

    initVector(Weight,-0.1,0.1);

    initVector(center,-4.0,4.0);

    initVector(delta,0.1,0.3);

    /*产生输入样本*/

    generateSample();

    /*开始迭代*/

    int iteration=ITERATION_CEIL;

    while(iteration-->0){

        if(calTotalError()<ERR)      //误差已达到要求，可以退出迭代

            break;

        updateParam();      //更新网络参数

    }

    cout<<"迭代次数:"<<ITERATION_CEIL-iteration-1<<endl;

     

    //根据已训练好的神经网络作几组测试

    for(int x=-4;x<5;++x){

        cout<<x<<" ";

        cout<<setprecision(8)<<setiosflags(ios::left)<<setw(15);

        cout<<getOutput(x)<<Hermit(x)<<endl;      //先输出我们预测的值，再输出真实值

    }

    return 0;

}

原文来自:博客园（华夏35度）http://www.cnblogs.com/zhangchaoyang 作者:Orisun

 

径向基网络（RBF network）之BP监督训练

径向基网络（RBF network）之BP监督训练

zouxy09@qq.com

http://blog.csdn.net/zouxy09

       之前看了流行学习的时候，感觉它很神奇，可以将一个4096维的人脸图像降到3维。然后又看到了可以用径向基网络来将这3维的图像重构到4096维。看到效果的时候，我和小伙伴们都惊呆了（呵呵，原谅我的孤陋寡闻）。见下图，第1和3行是原图像，维度是64x64=4096维，第2和第4行是将4096维的原图像用流行学习算法降到3维后，再用RBF网络重构回来的图像（代码是参考一篇论文写的）。虽然在重构领域，这效果不一定是好的，但对于无知的我，其中的奥妙勾引了我，使我忍不住又去瞻仰了一番。

       推荐大家先看看这个博主的这篇博文：

http://www.cnblogs.com/zhangchaoyang/articles/2591663.html

一、径向基函数

       在说径向基网络之前，先聊下径向基函数（Radical Basis Function，RBF）。径向基函数（Radical Basis Function，RBF）方法是Powell在1985年提出的。所谓径向基函数，其实就是某种沿径向对称的标量函数。通常定义为空间中任一点x到某一中心c之间欧氏距离的单调函数，可记作k(||x-c||)，其作用往往是局部的，即当x远离c时函数取值很小。例如高斯径向基函数：

       当年径向基函数的诞生主要是为了解决多变量插值的问题。可以看下面的图。具体的话是先在每个样本上面放一个基函数，图中每个蓝色的点是一个样本，然后中间那个图中绿色虚线对应的，就表示的是每个训练样本对应一个高斯函数（高斯函数中心就是样本点）。然后假设真实的拟合这些训练数据的曲线是蓝色的那根（最右边的图），如果我们有一个新的数据x1，我们想知道它对应的f(x1)是多少，也就是a点的纵坐标是多少。那么由图可以看到，a点的纵坐标等于b点的纵坐标加上c点的纵坐标。而b的纵坐标是第一个样本点的高斯函数的值乘以一个大点权值得到的，c的纵坐标是第二个样本点的高斯函数的值乘以另一个小点的权值得到。而其他样本点的权值全是0，因为我们要插值的点x1在第一和第二个样本点之间，远离其他的样本点，那么插值影响最大的就是离得近的点，离的远的就没什么贡献了。所以x1点的函数值由附近的b和c两个点就可以确定了。拓展到任意的新的x，这些红色的高斯函数乘以一个权值后再在对应的x地方加起来，就可以完美的拟合真实的函数曲线了。

二、径向基网络

       到了1988年， Moody和 Darken提出了一种神经网络结构，即RBF神经网络，属于前向神经网络类型，它能够以任意精度逼近任意连续函数，特别适合于解决分类问题。

       RBF网络的结构与多层前向网络类似，它是一种三层前向网络。输入层由信号源结点组成；第二层为隐含层，隐单元数视所描述问题的需要而定，隐单元的变换函数是RBF径向基函数，它是对中心点径向对称且衰减的非负非线性函数；第三层为输出层，它对输入模式的作用作出响应。从输人空间到隐含层空间的变换是非线性的，而从隐含层空间到输出层空间变换是线性的。

       RBF网络的基本思想是：用RBF作为隐单元的“基”构成隐含层空间，这样就可将输入矢量直接（即不需要通过权连接）映射到隐空间。根据Cover定理，低维空间不可分的数据到了高维空间会更有可能变得可分。换句话来说，RBF网络的隐层的功能就是将低维空间的输入通过非线性函数映射到一个高维空间。然后再在这个高维空间进行曲线的拟合。它等价于在一个隐含的高维空间寻找一个能最佳拟合训练数据的表面。这点与普通的多层感知机MLP是不同的。

       当RBF的中心点确定以后，这种映射关系也就确定了。而隐含层空间到输出空间的映射是线性的，即网络的输出是隐单元输出的线性加权和，此处的权即为网络可调参数。由此可见，从总体上看，网络由输人到输出的映射是非线性的，而网络输出对可调参数而言却又是线性的。这样网络的权就可由线性方程组直接解出，从而大大加快学习速度并避免局部极小问题。

       从另一个方面也可以这样理解，多层感知器（包括BP神经网络）的隐节点基函数采用线性函数，激活函数则采用Sigmoid函数或硬极限函数。而RBF网络的隐节点的基函数采用距离函数（如欧氏距离），并使用径向基函数（如Gaussian函数）作为激活函数。径向基函数关于n维空间的一个中心点具有径向对称性，而且神经元的输入离该中心点越远，神经元的激活程度就越低。隐节点的这一特性常被称为“局部特性”。

三、RBF网络的设计与求解

       RBF的设计主要包括两个方面，一个是结构设计，也就是说隐藏层含有几个节点合适。另一个就是参数设计，也就是对网络各参数进行求解。由上面的输入到输出的网络映射函数公式可以看到，网络的参数主要包括三种：径向基函数的中心、方差和隐含层到输出层的权值。到目前为止，出现了很多求解这三种参数的方法，主要可以分为以下两大类：

1、方法一：

       通过非监督方法得到径向基函数的中心和方差，通过监督方法（最小均方误差）得到隐含层到输出层的权值。具体如下：

（1）在训练样本集中随机选择h个样本作为h个径向基函数的中心。更好的方法是通过聚类，例如K-means聚类得到h个聚类中心，将这些聚类中心当成径向基函数的h个中心。

（2）RBF神经网络的基函数为高斯函数时，方差可由下式求解：

       式中c_max 为所选取中心之间的最大距离，h是隐层节点的个数。扩展常数这么计算是为了避免径向基函数太尖或太平。

（3）隐含层至输出层之间神经元的连接权值可以用最小均方误差LMS直接计算得到，计算公式如下：（计算伪逆）（d是我们期待的输出值）

2、方法二：

         采用监督学习算法对网络所有的参数（径向基函数的中心、方差和隐含层到输出层的权值）进行训练。主要是对代价函数（均方误差）进行梯度下降，然后修正每个参数。具体如下：

（1）随机初始化径向基函数的中心、方差和隐含层到输出层的权值。当然了，也可以选用方法一中的（1）来初始化径向基函数的中心。

（2）通过梯度下降来对网络中的三种参数都进行监督训练优化。代价函数是网络输出和期望输出的均方误差：

然后每次迭代，在误差梯度的负方向已一定的学习率调整参数。

四、代码实现：

1、第一种方法

       第一种方法在zhangchaoyang的博客上面有C++的实现，只是上面针对的是标量的数据（输入和输出都是一维的）。而在Matlab中也提供了第一种方法的改进版（呵呵，个人觉得，大家可以在Matlab中运行open newrb查看下源代码）。

       Matlab提供的一个函数是newrb()。它有个技能就是可以自动增加网络的隐层神经元数目直到均方差满足我们要求的精度或者神经元数数目达到最大（也就是我们提供的样本数目，当神经元数目和我们的样本数目一致时，rbf网络此时的均方误差为0）为止。它使用方法也能简单：

rbf = newrb(train_x, train_y);

output = rbf(test_x);

      直接把训练样本给它就可以得到一个rbf网络了。然后我们把输入给它就可以得到网络的输出了。

2、第二种方法

       第二种方法在zhangchaoyang的博客上面也有C++的实现，只是上面针对的还是标量的数据（输入和输出都是一维的）。但我是做图像的，网络需要接受高维的输入，而且在Matlab中，向量的运算要比for训练的运算要快很多。所以我就自己写了个可以接受向量输入和向量输出的通过BP算法监督训练的版本。BP算法可以参考这里：BackpropagationAlgorithm ，主要是计算每层每个节点的残差就可以了。另外，我的代码是可以通过梯度检查的，但在某些训练集上面，代价函数值却会随着迭代次数上升，这就很奇怪了，然后降低了学习率还是一样。但在某些简单点的训练集上面还是可以工作的，虽然训练误差也挺大的（没有完全拟合训练样本）。所以大家如果发现代码里面有错误的部分，还望大家告知下。

       主要代码见下面：

learnRBF.m

%// This is a RBF network trained by BP algorithm
%// Author : zouxy
%// Date   : 2013-10-28
%// HomePage : http://blog.csdn.net/zouxy09
%// Email  : zouxy09@qq.com

close all; clear; clc;

%%% ************************************************
%%% ************ step 0: load data ****************
display('step 0: load data...');
% train_x = [1 2 3 4 5 6 7 8]; % each sample arranged as a column of train_x
% train_y = 2 * train_x;
train_x = rand(5, 10);
train_y = 2 * train_x;
test_x = train_x;
test_y = train_y;

%% from matlab
% rbf = newrb(train_x, train_y);
% output = rbf(test_x);


%%% ************************************************
%%% ******** step 1: initialize parameters ********
display('step 1: initialize parameters...');
numSamples = size(train_x, 2);
rbf.inputSize = size(train_x, 1);
rbf.hiddenSize = numSamples;        % num of Radial Basis function
rbf.outputSize = size(train_y, 1);
rbf.alpha = 0.1;  % learning rate (should not be large!)

%% centre of RBF
for i = 1 : rbf.hiddenSize
    % randomly pick up some samples to initialize centres of RBF
    index = randi([1, numSamples]);
    rbf.center(:, i) =  train_x(:, index);
end

%% delta of RBF
rbf.delta = rand(1, rbf.hiddenSize);

%% weight of RBF
r = 1.0; % random number between [-r, r]
rbf.weight = rand(rbf.outputSize, rbf.hiddenSize) * 2 * r - r;


%%% ************************************************
%%% ************ step 2: start training ************
display('step 2: start training...');
maxIter = 400;
preCost = 0;
for i = 1 : maxIter
    fprintf(1, 'Iteration %d ,', i);
    rbf = trainRBF(rbf, train_x, train_y);
    fprintf(1, 'the cost is %d 
', rbf.cost);

    curCost = rbf.cost;
    if abs(curCost - preCost) < 1e-8
        disp('Reached iteration termination condition and Termination now!');
        break;
    end
    preCost = curCost;
end


%%% ************************************************
%%% ************ step 3: start testing ************
display('step 3: start testing...');
Green = zeros(rbf.hiddenSize, 1);
for i = 1 : size(test_x, 2)
    for j = 1 : rbf.hiddenSize
        Green(j, 1) = green(test_x(:, i), rbf.center(:, j), rbf.delta(j));
    end
    output(:, i) = rbf.weight * Green;
end
disp(test_y);
disp(output);

trainRBF.m

function [rbf] = trainRBF(rbf, train_x, train_y)
    %%% step 1: calculate gradient
    numSamples = size(train_x, 2);
    Green = zeros(rbf.hiddenSize, 1);
    output = zeros(rbf.outputSize, 1);
    delta_weight = zeros(rbf.outputSize, rbf.hiddenSize);
    delta_center = zeros(rbf.inputSize, rbf.hiddenSize);
    delta_delta =  zeros(1, rbf.hiddenSize);
    rbf.cost = 0;
    for i = 1 : numSamples
        %% Feed forward
        for j = 1 : rbf.hiddenSize
            Green(j, 1) = green(train_x(:, i), rbf.center(:, j), rbf.delta(j));
        end
        output = rbf.weight * Green;

        %% Back propagation
        delta3 = -(train_y(:, i) - output);
        rbf.cost = rbf.cost + sum(delta3.^2);
        delta_weight = delta_weight + delta3 * Green';
        delta2 = rbf.weight' * delta3 .* Green;
        for j = 1 : rbf.hiddenSize
            delta_center(:, j) = delta_center(:, j) + delta2(j) .* (train_x(:, i) - rbf.center(:, j)) ./ rbf.delta(j)^2;
            delta_delta(j) = delta_delta(j)+ delta2(j) * sum((train_x(:, i) - rbf.center(:, j)).^2) ./ rbf.delta(j)^3;
        end
    end

    %%% step 2: update parameters
    rbf.cost = 0.5 * rbf.cost ./ numSamples;
    rbf.weight = rbf.weight - rbf.alpha .* delta_weight ./ numSamples;
    rbf.center = rbf.center - rbf.alpha .* delta_center ./ numSamples;
    rbf.delta = rbf.delta - rbf.alpha .* delta_delta ./ numSamples;
end

green.m

function greenValue = green(x, c, delta)
    greenValue = exp(-1.0 * sum((x - c).^2) / (2 * delta^2));
end

五、代码测试

      首先，我测试了一维的输入，需要拟合的函数很简单，就是y=2x。

train_x = [1 2 3 4 5 6 7 8];

train_y = 2 * train_x;

     所以期待的输出就是：

2    4     6     8   10    12    14   16

     我代码训练迭代200次后的网络输出是：

2.0042   4.0239    5.9250    8.0214  10.0692   11.9351   14.0179  15.9958

      Matlab的newrb的输出是：

2.0000   4.0000    6.0000    8.0000  10.0000   12.0000   14.0000  16.0000

      可以看到，Matlab的是完美拟合啊。我的那个还是均方误差还是挺大的。

      然后，我测试了高维的输入，训练样本是通过Matlab的rand(5, 10)来得到的，它生成的是5行10列[0 1]之间的随机数。也就是说我们的样本是10个，每个样本的维度是5维。我们测试的也是很简单的函数y=2x。结果如下：

 

      关于这个结果，我也不说什么了。期待大家发现代码里面错误的地方，然后告知下，非常感谢。

RBF神经网络与BP神经网络的比较

RBF神经网络与BP神经网络的比较

RBF神经网络与BP神经网络都是非线性多层前向网络，它们都是通用逼近器。对于任一个BP神经网络，总存在一个RBF神经网络可以代替它，反之亦然。但是这两个网络也存在着很多不同点，这里从网络结构、训练算法、网络资源的利用及逼近性能等方面对RBF神经网络和BP神经网络进行比较研究。

（1） 从网络结构上看。 BP神经网络实行权连接，而RBF神经网络输入层到隐层单元之间为直接连接，隐层到输出层实行权连接。BP神经网络隐层单元的转移函数一般选择非线性函数（如反正切函数），RBF神经网络隐层单元的转移函数是关于中心对称的RBF（如高斯函数）。BP神经网络是三层或三层以上的静态前馈神经网络，其隐层和隐层节点数不容易确定，没有普遍适用的规律可循，一旦网络的结构确定下来，在训练阶段网络结构将不再变化；RBF神经网络是三层静态前馈神经网络，隐层单元数也就是网络的结构可以根据研究的具体问题，在训练阶段自适应地调整，这样网络的适用性就更好了。

（2） 从训练算法上看。 BP神经网络需要确定的参数是连接权值和阈值，主要的训练算法为BP算法和改进的BP算法。但BP算法存在许多不足之处，主要表现为易限于局部极小值，学习过程收敛速度慢，隐层和隐层节点数难以确定；更为重要的是，一个新的BP神经网络能否经过训练达到收敛还与训练样本的容量、选择的算法及事先确定的网络结构（输入节点、隐层节点、输出节点及输出节点的传递函数）、期望误差和训练步数有很大的关系。RBF神经网络的训练算法在前面已做了论述，目前，很多RBF神经网络的训练算法支持在线和离线训练，可以动态确定网络结构和隐层单元的数据中心和扩展常数，学习速度快，比BP算法表现出更好的性能。

（3） 从网络资源的利用上看。 RBF神经网络原理、结构和学习算法的特殊性决定了其隐层单元的分配可以根据训练样本的容量、类别和分布来决定。如采用最近邻聚类方式训练网络，网络隐层单元的分配就仅与训练样本的分布及隐层单元的宽度有关，与执行的任务无关。在隐层单元分配的基础上，输入与输出之间的映射关系，通过调整隐层单元和输出单元之间的权值来实现，这样，不同的任务之间的影响就比较小，网络的资源就可以得到充分的利用。这一点和BP神经网络完全不同，BP神经网络权值和阈值的确定由每个任务（输出节点）均方差的总和直接决定，这样，训练的网络只能是不同任务的折中，对于某个任务来说，就无法达到最佳的效果。而RBF神经网络则可以使每个任务之间的影响降到较低的水平，从而每个任务都能达到较好的效果，这种并行的多任务系统会使RBF神经网络的应用越来越广泛。

总之，RBF神经网络可以根据具体问题确定相应的网络拓扑结构，具有自学习、自组织、自适应功能，它对非线性连续函数具有一致逼近性，学习速度快，可以进行大范围的数据融合，可以并行高速地处理数据。RBF神经网络的优良特性使得其显示出比BP神经网络更强的生命力，正在越来越多的领域内替代BP神经网络。目前，RBF神经网络已经成功地用于非线性函数逼近、时间序列分析、数据分类、模式识别、信息处理、图像处理、系统建模、控制和故障诊断等。

rbf原理 

所谓径向基函数 (Radial Basis Function 简称 RBF), 就是某种沿径向对称的标量函数。 
通常定义为空间中任一点x到某一中心xc之间欧氏距离的单调函数 , 可记作 k(||x-xc||), 
其作用往往是局部的 , 即当x远离xc时函数取值很小。最常用的径向基函数是高斯核函数 ,
形式为 k(||x-xc||)=exp{- ||x-xc||^2/(2*σ)^2) } 其中xc为核函数中心,σ为函数的宽度参数 , 
控制了函数的径向作用范围。在RBF网络中,这两个参数往往是可调的。 

可以从两个方面理解 RBF 网络的作用 : 
(1)把网络看成对未知函数f(x)的逼近器。
   一般任何函数都可表示成一组基函数的加权和 ,这相当于用隐层单元的输出函数构成一组基函数来逼近f(x) 
(2)在RBF网络中以输入层到隐层的基函数输出是一种非线性映射，而输出则是线性的。
   这样,RBF网络可以看成是首先将原始的非线性可分的特征空间变换到另一空间（通常是高维空间）,
   通过合理选择这一变换使在新空间中原问题线性可分,然后用一个线性单元元来解决问题。 
在典型的RBE网络中有三组可调参数：隐层基函数中心、方差,以及输出单元的权值。

这些参数的选择有三种常见的方法：
 (1)根据经验选择函数中心。
   比如只要训练样本的分布能代表所给问题 ,可根据经验选定均匀分布的M个中心, 
   其间距为d，可选取高斯核函数的方为σ=d/sqrt(2*M)。
 (2)用聚类方法选择基函数。
    可以各聚类中心作为核函数中心,而以各类样本的方差的某一函数作为各个基函数的宽度参数。 
   用(1)或(2)的方法选定了隐层基函旗的参数后,因输出单元是线性单元,它的权值可以简单地用最小二乘法
   直接计算出来。
 (3)将三组可调参数都通过训练样本用误差纠正算法求得。
   做法与BP方法类似,分别计算误差e(k)对各组参数的偏导数，然后用迭代求取参数。 
   研究表明,用于模式识别问题的RBF网络在一定意义上等价于首先用非参数方法估计出概率密度，
   必然后用它进行分类
http://www.2nsoft.cn/bbs/read.php?tid=741&fpage=2

RBF神经网络