极限学习机(ELM)代码注解

此代码来自黄老师主页，测试数据也可在主页下载。另外，代码形式略微改动，但本质未变！原理可查看：ELM原理。

---

function [TrainingTime, TestingTime, TrainingAccuracy, TestingAccuracy] = elm(TrainingData_File, TestingData_File, Elm_Type, NumberofHiddenNeurons, ActivationFunction)

% Usage: elm(TrainingData_File, TestingData_File, Elm_Type, NumberofHiddenNeurons, ActivationFunction)
% OR:    [TrainingTime, TestingTime, TrainingAccuracy, TestingAccuracy] = elm(TrainingData_File, TestingData_File, Elm_Type, NumberofHiddenNeurons, ActivationFunction)
%
% Input:
% TrainingData_File     - Filename of training data set
% TestingData_File      - Filename of testing data set
% Elm_Type              - 0 for regression; 1 for (both binary and multi-classes) classification
% NumberofHiddenNeurons - Number of hidden neurons assigned to the ELM
% ActivationFunction    - Type of activation function:
%                           'sig' for Sigmoidal function
%                           'sin' for Sine function
%                           'hardlim' for Hardlim function
%                           'tribas' for Triangular basis function
%                           'radbas' for Radial basis function (for additive type of SLFNs instead of RBF type of SLFNs)
%
% Output: 
% TrainingTime          - Time (seconds) spent on training ELM
% TestingTime           - Time (seconds) spent on predicting ALL testing data
% TrainingAccuracy      - Training accuracy: 
%                           RMSE for regression or correct classification rate for classification
% TestingAccuracy       - Testing accuracy: 
%                           RMSE for regression or correct classification rate for classification
%
% MULTI-CLASSE CLASSIFICATION: NUMBER OF OUTPUT NEURONS WILL BE AUTOMATICALLY SET EQUAL TO NUMBER OF CLASSES
% FOR EXAMPLE, if there are 7 classes in all, there will have 7 output
% neurons; neuron 5 has the highest output means input belongs to 5-th class
%
% Sample1 regression: [TrainingTime, TestingTime, TrainingAccuracy, TestingAccuracy] = elm('sinc_train', 'sinc_test', 0, 20, 'sig')
% Sample2 classification: elm('diabetes_train', 'diabetes_test', 1, 20, 'sig')
%
    %%%%    Authors:    MR QIN-YU ZHU AND DR GUANG-BIN HUANG
    %%%%    NANYANG TECHNOLOGICAL UNIVERSITY, SINGAPORE
    %%%%    EMAIL:      EGBHUANG@NTU.EDU.SG; GBHUANG@IEEE.ORG
    %%%%    WEBSITE:    http://www.ntu.edu.sg/eee/icis/cv/egbhuang.htm
    %%%%    DATE:       APRIL 2004
    
%%%-----------------------------本程序（函数）算法摘要：----------------------------------------------------------%%%
    %训练样本数：N  测试样本数：N'  输入神经元：n  输出神经元：m  隐层神经元：L
    %输入层→隐层：    H=g(W·X+b)   (L×N=L×n×n×N +L×N)
    %训练网络参数β：  β=(H')+·T’  （L×m=L×N×N×m）
    %计算实际输出：    Y=(H'·β)'   （m×N=(N×L×L×m)'）
    %回归时：m=1(仅1类输出神经元,即回归值为行向量)   分类时：m=m（m即输出神经元数，来自训练、测试总和label中）
    %如果是回归：训练、测试精度越小越好！   如果是分类：训练、测试准确度越大越好！

%%%-------------------------宏定义：elm的类型---------------------------------------------------------%%%
%%%%%%%%%%% Macro definition
REGRESSION=0;
CLASSIFIER=1;

%%%-----------------------------------数据加载-------------------------------------------------------%%%
%%%%%%%%%%% Load training dataset
train_data=load(TrainingData_File);
T=train_data(:,1)';    %（期望输出→回归训练无需处理，分类训练需处理） T：“第一列”并转置（1×N）
P=train_data(:,2:size(train_data,2))';    %（训练数据）P：剩下所有的列并转置（n×N,n为输入神经元个数）
clear train_data;                               %   Release raw training data array

%%%%%%%%%%% Load testing dataset
test_data=load(TestingData_File);
TV.T=test_data(:,1)';    %（期望输出→回归测试无需处理，分类测试需处理） TV.T：第一列并转置(1×N')
TV.P=test_data(:,2:size(test_data,2))';   %（测试数据）TV.P：剩下的列并转置（n'×N'=n×N'）
clear test_data;                                    %   Release raw testing data array

NumberofTrainingData=size(P,2); %训练样本数（N）
NumberofTestingData=size(TV.P,2); %测试数据组数(N')
NumberofInputNeurons=size(P,1); %输入神经元个数(n,训练=测试)

%如前所述：X(P)、T已经构造完毕；相关参数：L、n(训练=测试,因同用一个网络)、N、N'已知
%如果作回归分析，则可直接计算（获得输出权重β(L×1),再计算实际输出Y(1×N)）
%如果需作分类,则先计算输出神经元m,再构造训练、测试各自的期望输出T
%%%----------------------------------------------------------------------------------------------%%%


%%%----------------------------分类预处理（仅用于分类）--------------------------------------------%%%%
if Elm_Type == CLASSIFIER 
    %%%%%%%%%%%% Preprocessing the data of classification
    sorted_target=sort(cat(2,T,TV.T),2);  %先水平连接训练T、测试TV.T（得1×(N+N')）,再按升序对其“行排列“
    label=zeros(1,1);                            %   Find and save in 'label' class label from training and testing data sets
    label(1,1)=sorted_target(1,1);
    j=1;
    for i = 2:(NumberofTrainingData+NumberofTestingData)  
        if sorted_target(1,i) ~= label(1,j)  %去除sorted_target（训练+测试）中重复的标签→存入label（1×m）
            j=j+1;      %label包含：训练、测试
            label(1,j) = sorted_target(1,i);
        end
    end
    number_class=j;
    NumberofOutputNeurons=number_class; 
       
    %%%%%%%%%% Processing the targets of training   构造分类的训练期望输出T(m×N)
    temp_T=zeros(NumberofOutputNeurons, NumberofTrainingData);  %temp_T(m×N)
    for i = 1:NumberofTrainingData    %外：列（横向→N）
        for j = 1:number_class     %内：行（竖向→m）
            if label(1,j) == T(1,i)  
                break; 
            end
        end
        temp_T(j,i)=1; %对temp_T每一列来说（T(1,i)对应temp_T第i列），标签T(1,i)匹配了label(1×m)中第j个标签,则把temp_T此列的此元素（即temp_T(j,i)）赋值为1（仅赋这第个）
    end
    T=temp_T*2-1;  %把temp_T中（1,0）→（1，-1）

    %%%%%%%%%% Processing the targets of testing   构造分类的测试期望输出T(m×N')
    temp_TV_T=zeros(NumberofOutputNeurons, NumberofTestingData);  %（m×N'）
    for i = 1:NumberofTestingData
        for j = 1:number_class
            if label(1,j) == TV.T(1,i)
                break; 
            end
        end
        temp_TV_T(j,i)=1;
    end
    TV.T=temp_TV_T*2-1;
end                                           %   end if of Elm_Type
%分类预处理仅仅使得：分类的T≠回归的T
%%%---------------------------分类预处理结束--------------------------------------------------%%%

%%%-------------------------先训练，后测试（回归与分类通用）------------------------------------%%%
%（1）训练网络
%%%%%%%%%%% Calculate weights & biases
start_time_train=cputime;  %开始计时

%%%%%%%%%%% Random generate input weights InputWeight (w_i) and biases BiasofHiddenNeurons (b_i) of hidden neurons
InputWeight=rand(NumberofHiddenNeurons,NumberofInputNeurons)*2-1;  %InputWeigth(L×n),[0,1]均匀分布扩展至[-1,1]
BiasofHiddenNeurons=rand(NumberofHiddenNeurons,1);  %b(L×1)
tempH=InputWeight*P;  %tempH=W·X (L×N=L×n×n×N)
clear P;                                            %   Release input of training data 
ind=ones(1,NumberofTrainingData);
BiasMatrix=BiasofHiddenNeurons(:,ind); %b(L×N)     %   Extend the bias matrix BiasofHiddenNeurons to match the demention of H
tempH=tempH+BiasMatrix; %tempH=W·X+b

%%%%%%%%%%% Calculate hidden neuron output matrix H
switch lower(ActivationFunction)  
    case {'sig','sigmoid'}
        %%%%%%%% Sigmoid 
        H = 1 ./ (1 + exp(-tempH));  %H(L×N)=g(temp_H)
    case {'sin','sine'}
        %%%%%%%% Sine
        H = sin(tempH);    
    case {'hardlim'}
        %%%%%%%% Hard Limit
        H = double(hardlim(tempH));
    case {'tribas'}
        %%%%%%%% Triangular basis function
        H = tribas(tempH);
    case {'radbas'}
        %%%%%%%% Radial basis function
        H = radbas(tempH);
        %%%%%%%% More activation functions can be added here                
end
clear tempH;                                        %   Release the temparary array for calculation of hidden neuron output matrix H

%%%%%%%%%%% Calculate output weights OutputWeight (beta_i)
OutputWeight=pinv(H') * T';   %OutputWeight(L×m=L×N × N×m)     %（无正则化因子）implementation without regularization factor //refer to 2006 Neurocomputing paper
%%%T(m×N)、H(L×N),程序欲使用（Hβ=T→β=H+·T）→则H、T都需要转置：H'*β=T'(N×L×L×m=N×m)！

end_time_train=cputime;  %结束计时
TrainingTime=end_time_train-start_time_train;       %（训练时间）Calculate CPU time (seconds) spent for training ELM

%%%%%%%%%%% Calculate the training accuracy
Y=(H' * OutputWeight)';  %m×N=(N×L×L×m)'          % (训练的实际输出)Y: the actual output of the training data
if Elm_Type == REGRESSION
    TrainingAccuracy=sqrt(mse(T - Y));  %回归训练均方根误差(越小越小！)     % Calculate training accuracy (RMSE) for regression case
elseif Elm_Type == CLASSIFIER
    MissClassificationRate_Training=0;
    for i = 1 : size(T, 2)  %T/Y(m×N)
        [x, label_index_expected]=max(T(:,i));
        [x, label_index_actual]=max(Y(:,i));
        if label_index_actual~=label_index_expected
            MissClassificationRate_Training=MissClassificationRate_Training+1;
        end
    end
    TrainingAccuracy=1-MissClassificationRate_Training/size(T,2);  %分类训练准确度（越大越好！）
end
clear H;

%（2）测试网络
%前面网络已训练完毕→获得β(OutputWeigth)；欲用此网络处理测试数据,仅需测试数据自己的X、b→H→TY
%%%%%%%%%%% Calculate the output of testing input
start_time_test=cputime;
tempH_test=InputWeight*TV.P;   %InputWeigth(L×n)、TV.P(n×N')
clear TV.P;            %   Release input of testing data             
ind=ones(1,NumberofTestingData); %(1×N')
BiasMatrix=BiasofHiddenNeurons(:,ind);   %(L×N')          %   Extend the bias matrix BiasofHiddenNeurons to match the demention of H
tempH_test=tempH_test + BiasMatrix;
switch lower(ActivationFunction)
    case {'sig','sigmoid'}
        %%%%%%%% Sigmoid 
        H_test = 1 ./ (1 + exp(-tempH_test));  %H_test(L×N') 
    case {'sin','sine'}
        %%%%%%%% Sine
        H_test = sin(tempH_test);        
    case {'hardlim'}
        %%%%%%%% Hard Limit
        H_test = hardlim(tempH_test);        
    case {'tribas'}
        %%%%%%%% Triangular basis function
        H_test = tribas(tempH_test);        
    case {'radbas'}
        %%%%%%%% Radial basis function
        H_test = radbas(tempH_test);        
        %%%%%%%% More activation functions can be added here        
end
TY=(H_test' * OutputWeight)';  %(m×N'=(N'×L×L×m)')    % TY: the actual output of the testing data
end_time_test=cputime;
TestingTime=end_time_test-start_time_test;           %   Calculate CPU time (seconds) spent by ELM predicting the whole testing data

if Elm_Type == REGRESSION
    TestingAccuracy=sqrt(mse(TV.T - TY));           %  Calculate testing accuracy (RMSE) for regression case
elseif Elm_Type ==CLASSIFIER
    MissClassificationRate_Testing=0;
    for i = 1 : size(TV.T, 2)   %TV.T/TY(m×N')
        [x, label_index_expected]=max(TV.T(:,i));
        [x, label_index_actual]=max(TY(:,i));
        if label_index_actual~=label_index_expected
            MissClassificationRate_Testing=MissClassificationRate_Testing+1;
        end
    end
    TestingAccuracy=1-MissClassificationRate_Testing/size(TV.T,2);  
end
end

---