徑向基函數（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<< "\t" ;

         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<< "\t" ;

         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 \n', 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神經網絡