徑向基函數（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[]){
　　 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;
}

（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;
}