本博客来源于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进行改进，更多内容请点击专栏获取。
)
