本博客來源于CSDN機器魚，未同意任何人轉載。
更多內容，歡迎點擊本專欄目錄，查看更多內容。
目錄
0.引言
1.原理
2.具體實現
3.結語
0.引言
在博客【深度核極限學習機】里我們講述了深度核極限學習機原理，今天我們對其繼續進行改進：觀察下列DKELM的調用代碼可以發現，超參數有隱含層節點數、ELM-AE的L2正則化系數，核函數的懲罰系數與核參數等，這些參數都直接影響分類或回歸的準確率，做過SVM的都知道，光懲罰系數與核參數的選擇都很困難，一般也不是自己手動選擇，都是用各種優化智能優化算法進行手動選擇，因此我們繼續利用優化算法進行優化選擇。
%elm-ae的參數
h=[100 ,50];%各隱含層節點數，里面有幾個數就是幾個隱含層，即有幾個ELM-AE參與無監督預訓練，數字代表隱含層節點數
% 比如 h=[100,50],假如輸入是200，輸出是5類，則DKELM網絡結構為200-100-50-5。
% 比如 h=[100,50，20],假如輸入是200，輸出是5類，則DKELM網絡結構為200-100-50-20-5。
TF='sig';%ELM-AE激活函數
lambda=1000;%ELM-AEL2正則化系數
%頂層核極限學習的懲罰系數與核參數
kernel='RBF_kernel';%核函數類型 RBF_kernel lin_kernel poly_kernel wav_kernel
% kernelpara中第一個參數是懲罰系數
if strcmp(kernel,'RBF_kernel')kernelpara=[1 1]; %rbf有一個核參數
elseif strcmp(kernel,'lin_kernel')kernelpara=[1]; %線性沒有核參數
elseif strcmp(kernel,'poly_kernel')kernelpara=[1,1,1]; %多項式有2個核參數
elseif strcmp(kernel,'wav_kernel')kernelpara=[1,1,1,1]; %小波核有3個核參數
end
%%
delm=dkelmtrain(P_train,T_train,h,lambda,TF,kernelpara,kernel);
T1=dkelmpredict(delm,kernelpara,kernel,P_train,P_train);1.原理
麻雀優化算法的原理到處都是，就不詳細說明了。我們只需要知道一點，任何一個優化算法，都是在他的尋優空間內，基于特定公式進行尋優，具體什么公式我們可以不用管，只要網上能找到可以跑得通的原始代碼，我們拿過來改一下就能用，我要調整的就是尋優參數的維度，每個參數的尋優范圍，以及適應度函數。
具體實現如下：
2.具體實現
步驟1：確定尋優參數與尋優的范圍。如下代碼所示，我們要優化的是各隱含層節點數，ELM-AE的正則化系數，頂層kelm的懲罰系數與核參數，隱含層數量與核函數類型是預先設定好的，一般設置成3層與RBF核函數即可。由于每個參數的代碼的含義不一樣，因此我們需要針對每個參數設定尋優范圍lb與ub，詳細的看下列代碼。
function [bestX,Convergence_curve,result]=ssafordkelm(P_train,P_test,T_train,T_test,layers,TF,kernel)
% 優化DKELM超參數,包括各隱含層節點數，正則化系數，頂層kelm的懲罰系數與核參數
% 定義尋優邊界，由于隱含層層數與核函數類型會影響尋優維度，所以下面我們分段進行尋優范圍設置
lb=[];ub=[];
% 首先是各層節點數,范圍是1-100，整數
lb=[lb ones(1,layers)];
ub=[ub 100*ones(1,layers)];
% 然后是ELM-AE的正則化系數與 kelm的懲罰系數
lb=[lb 1e-3 1e-3];
ub=[ub 1e3 1e3];
% 最后是kelm的核參數
% 核參數設置 詳情看kernel_matrix
if strcmp(kernel,'RBF_kernel')lb=[lb 1e-3];ub=[ub 1e3 ]; %rbf有一個核參數
elseif strcmp(kernel,'lin_kernel')lb=[lb];ub=[ub]; %線性沒有核參數
elseif strcmp(kernel,'poly_kernel')lb=[lb 1e-3 1];ub=[ub 1e3 10]; %多項式有2個核參數,第二個是多項式冪指數，我們把范圍就定到1-10
elseif strcmp(kernel,'wav_kernel')lb=[lb 1e-3 1e-3 1e-3];ub=[ub 1e3 1e3 1e3]; %小波核有3個核參數
end
dim=length(lb);pop=10; % 種群數
M=10; % Maximum numbef of iterationsP_percent = 0.2; % The population size of producers accounts for "P_percent" percent of the total population size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pNum = round( pop * P_percent ); % The population size of the producers%Initialization
for i=1:pop%隨機初始化位置for j=1:dimif j>layers% % 隱含層節點是整數 其他的都是浮點數x(i,j)=(ub(j)-lb(j))*rand+lb(j);elsex(i,j)=round((ub(j)-lb(j))*rand+lb(j)); %endendfit(i)=fitness(x(i,:),P_train,T_train,P_test,T_test,layers,TF,kernel);
endpFit = fit;
pX = x;
fMin=fit(1);
bestX = x( i, : );
result=[];
for t = 1 : Mt[ ans, sortIndex ] = sort( pFit );% Sort.從小到大[fmax,B]=max( pFit );worse= x(B,:);r2=rand(1);%%%%%%%%%%%%%5%%%%%%這一部位為發現者（探索者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%if(r2<0.8)%預警值較小，說明沒有捕食者出現for i = 1 : pNum %r2小于0.8的發現者的改變（1-20） % Equation (3)r1=rand(1);x( sortIndex( i ), : ) = pX( sortIndex( i ), : )*exp(-(i)/(r1*M));%對自變量做一個隨機變換x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub,layers );%對超過邊界的變量進行去除fit( sortIndex( i ) )=fitness(x( sortIndex( i ), : ),P_train,T_train,P_test,T_test,layers,TF,kernel);endelse %預警值較大，說明有捕食者出現威脅到了種群的安全，需要去其它地方覓食for i = 1 : pNum %r2大于0.8的發現者的改變x( sortIndex( i ), : ) = pX( sortIndex( i ), : )+randn(1)*ones(1,dim);x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub,layers );fit( sortIndex( i ) )= fitness(x( sortIndex( i ), : ),P_train,T_train,P_test,T_test,layers,TF,kernel);endend[ fMMin, bestII ] = min( fit );bestXX = x( bestII, : );%%%%%%%%%%%%%5%%%%%%這一部位為加入者（追隨者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%for i = ( pNum + 1 ) : pop %剩下20-100的個體的變換 % Equation (4)
% i
% sortIndex( i )A=floor(rand(1,dim)*2)*2-1;if( i>(pop/2))%這個代表這部分麻雀處于十分饑餓的狀態（因為它們的能量很低，也是是適應度值很差），需要到其它地方覓食x( sortIndex(i ), : )=randn(1,dim).*exp((worse-pX( sortIndex( i ), : ))/(i)^2);else%這一部分追隨者是圍繞最好的發現者周圍進行覓食，其間也有可能發生食物的爭奪，使其自己變成生產者x( sortIndex( i ), : )=bestXX+(abs(( pX( sortIndex( i ), : )-bestXX)))*(A'*(A*A')^(-1))*ones(1,dim);endx( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub,layers );%判斷邊界是否超出fit( sortIndex( i ) )=fitness(x( sortIndex( i ), : ),P_train,T_train,P_test,T_test,layers,TF,kernel);end%%%%%%%%%%%%%5%%%%%%這一部位為意識到危險（注意這里只是意識到了危險，不代表出現了真正的捕食者）的麻雀的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%c=randperm(numel(sortIndex));%%%%%%%%%這個的作用是在種群中隨機產生其位置（也就是這部分的麻雀位置一開始是隨機的，意識到危險了要進行位置移動，%處于種群外圍的麻雀向安全區域靠攏，處在種群中心的麻雀則隨機行走以靠近別的麻雀）b=sortIndex(c(1:3));for j = 1 : length(b) % Equation (5)if( pFit( sortIndex( b(j) ) )>(fMin) ) %處于種群外圍的麻雀的位置改變x( sortIndex( b(j) ), : )=bestX+(randn(1,dim)).*(abs(( pX( sortIndex( b(j) ), : ) -bestX)));else %處于種群中心的麻雀的位置改變x( sortIndex( b(j) ), : ) =pX( sortIndex( b(j) ), : )+(2*rand(1)-1)*(abs(pX( sortIndex( b(j) ), : )-worse))/ ( pFit( sortIndex( b(j) ) )-fmax+1e-50);endx( sortIndex(b(j) ), : ) = Bounds( x( sortIndex(b(j) ), : ), lb, ub,layers );fit( sortIndex( b(j) ) )=fitness(x( sortIndex( b(j) ), : ),P_train,T_train,P_test,T_test,layers,TF,kernel);endfor i = 1 : popif ( fit( i ) < pFit( i ) )pFit( i ) = fit( i );pX( i, : ) = x( i, : );endif( pFit( i ) < fMin )fMin= pFit( i );bestX = pX( i, : );endendConvergence_curve(t)=fMin;result(t,:)=bestX;
end步驟2：確定適應度函數。觀察ssafordkelm函數中更新fmin與bestX可知，這個優化算法是找適應度函數最小時對應的x值，是一個最小化問題，因此我們的適應度函數也要寫成最小化誤差的方式。代碼如下，首先將麻雀的x傳進來，進行解碼，得到各隱含層節點數、L2正則化系數、頂層kelm的懲罰系數與核參數，然后用這些參數訓練、預測DELM，得到預測值，分類任務就計算分類錯誤率，回歸任務就計算均方差。
function accuracy=fitness(xx,P_train,T_train,P_test,T_test,layers,activatation,kernel)
rng('default')hidden=xx(1:layers); %各隱含層節點數
lambda=xx(layers+1); %L2正則化系數
para=xx(layers+2:end); %頂層kelm懲罰系數與核參數delm=dkelmtrain(P_train,T_train,hidden,lambda,activatation,para,kernel);
T2=dkelmpredict(delm,para,kernel,P_test,P_train);
% 如果是分類任務
[~ ,J2]=max(T2,[],2);
[~ ,J3]=max(T_test,[],2);
accuracy=1-sum(J2==J3)/length(J2); %以最小化分類錯誤率作為適應度值，優化算法的目的就是找到一組最優超參數，用這組超參數訓練 的網絡，可以使得使得網絡的錯誤率最低% 如果是回歸任務
%[m,m]=size(T2);
%T2=reshape(T2,[1,m*n]);
%T_test=reshape(T_test,[1,m*n]);%accuracy=mse(T2,T_test); %以最小化預測值與實際值的均方差作為適應度值，優化算法的目的就是找到一組最優超參數，用這組超參數訓練 的網絡，可以使得使得網絡的預測誤差最低rng(sum(100*clock))步驟3：在ssa將x傳進適應度函數之前，我們還需要做一下參數范圍的限制，目的1是怕x中的值不在lb、ub中，目的2是部分參數是整數，部分是小數，我們需要將x的值轉成對應的類型。
function x = Bounds( x, xmin, xmax,layers)
D=length(x);
for j=1:Dif j <=layersx(j)=round(x(j));%隱含層節點數是整數endif x(j)>xmax(j) | x(j)<xmin(j)if j>layersx(j)=(xmax(j)-xmin(j))*rand+xmin(j); %elsex(j)=round((xmax(j)-xmin(j))*rand+xmin(j)); %endend
end步驟4：主程序調用
%% 優化DKELM超參數,包括各隱含層節點數，正則化系數，頂層kelm的懲罰系數與核參數,
close all;clear;format compact;format short;clc;warning off
%%
load data
P_train = double(train_x) ;%30*1000
P_test = double(test_x);%30*100
T_train = double(train_Y);%5*1000
T_test = double(test_Y);%5*100
%% 預先設定幾個常用參數，這些設置的參數不參與尋優
layers=3;%layers個隱含層
TF='sig';%ELM-AE激活函數
%頂層核極限學習機參數
kernel='RBF_kernel';%核函數類型 RBF_kernel lin_kernel poly_kernel wav_kernel%%
[x,trace,result]=ssafordkelm(P_train,P_test,T_train,T_test,layers,TF,kernel); %優化DKELM超參數,包括各隱含層節點數，正則化系數，頂層kelm的懲罰系數與核參數figure
plot(trace)
title('適應度曲線')
xlabel('優化代數')
ylabel('適應度值')%% 將優化得到的參數解碼出來
hidden=x(1:layers) %各隱含層節點數
lambda=x(layers+1) %L2正則化系數
kernelpara=x(layers+2:end) %頂層kelm懲罰系數與核參數%% 重新訓練DKELM，看結果
rng('default')
delm=dkelmtrain(P_train,T_train,hidden,lambda,TF,kernelpara,kernel);
T1=dkelmpredict(delm,kernelpara,kernel,P_train,P_train);%訓練集預測結果
T2=dkelmpredict(delm,kernelpara,kernel,P_test,P_train);%測試集預測結果
%如果是分類就計算分類正確率，如果是回歸就計算預測的誤差
其中用到的dkelmtrain dkelmpredict等函數，看我上一篇博客。
3.結語
本章我們將DKELM與麻雀優化算法結合起來。用麻雀優化算法來對dklem的參數進行尋優，下一篇我們繼續對DKELM進行改進，更多內容請點擊專欄獲取。
)
