kernel lda.m

用MATLAB编写的人脸检测算法

clear all;
close all;
ratio = 2;
class_num = 40;
sample_num = 10;
row_0 = 112;
line_0 = 92;
row = ceil(row_0/ratio);
line = ceil(line_0 / ratio);
row_line = min(row, line);
rowline = row * line;
load orl_2;

disp('orl Two-phase KFD algorithm now begin');
right_num = [];
rr_0 = [];

for count = 2:6
    tic;
    count
    train_num = count;%每类样本训练数
    train_total = train_num * class_num;%40类样本总的训练数M
    test_num = sample_num - train_num;%每类样本测试数=样本数-训练数
    test_total = test_num * class_num;%40类样本总的测试数
   

    %%---step1:Perform KPCA in input space---%%
                  %%---求K---%%
    l_m=(ones(train_total,train_total))/train_total;
    
    orl_train = zeros(rowline, train_total);
    for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
            orl_temp = double(orl(:,:,j,i));
            for p = 1:row
                for q = 1:line
                    site = (p-1)*line + q;  
                    orl_train(site, kk) = orl_temp(p,q);
                end
            end
        end
    end 

    orl_test = zeros(rowline, test_total);
  	for i = 1:class_num
        for j = 1:test_num
            kk = (i-1)*test_num + j;
            orl_temp = double(orl(:,:,j+train_num,i));
            for p = 1:row
                for q = 1:line
                    site = (p-1)*line + q;  
                    orl_test(site, kk) = orl_temp(p,q);
                end
            end
        end
    end 

    K1=zeros(train_total,train_total);%利用orl_train求K~=k(Xi,Xj),Xi为训练样本中的特征向量
    K2=zeros(test_total,test_total);%利用orl_test求K~=k(Xi,Xj),Xi为测试样本中的特征向量
    for i=1:train_total
        for j=i:train_total
            %K1(i,j)=((orl_train(:,i))'*(orl_train(:,j))+1).^2;%选择多项式函数为kernel函数 
            K1(i,j)=(orl_train(:,i))'*(orl_train(:,j))+1;%选择线性函数为kernel函数 
            K1(j,i)=K1(i,j);
        end
    end
    for i=1:test_total
        for j=i:test_total
            %K2(i,j)=((orl_test(:,i))'*(orl_test(:,j))+1).^2;
            K2(i,j)=(orl_test(:,i))'*(orl_test(:,j))+1;
            K2(j,i)=K2(i,j);
        end
    end
    
    K=K1-l_m*K1-K1*l_m+l_m*K1*l_m;
    K_test=K2-l_m*K2-K2*l_m+l_m*K2*l_m;

    %%---求K的前m个非负特征值及其对应的特征向量，m为K的秩---%%
    m=rank(K) - class_num;
    
    [V_K, D_K] = eig(K);%求K的前m个大于0的最大本征值及其对应的标准正交本征向量
    %V_K1 = orth(V_K);
    dd = abs(eig(K));
    [dd_value, dd_site] = sort(dd);
    temp2 = m;
    temp1 = train_total - temp2 + 1;
    yy = [];
    Dt = zeros(temp2);
    jj = 0;
    for i = train_total:-1:temp1
        yy = [yy; V_K(:,dd_site(i))'];
        jj = jj + 1;
        Dt(jj,jj) = dd_value(i);   
    end
    yy = orth(yy');
    %yy = yy';
    D_K1=diag(Dt);

    %%利用等式(16)进行KPCA转换%%
    y_temp=zeros(train_total,temp2);
    for i=1:train_total
        for j=1:temp2
            y_temp(i,j)=yy(i,j)/sqrt(D_K1(j));
        end
    end
    for i=1:train_total
        y_train(:,i)=y_temp'*K(:,i);
        y_test(:,i)=y_temp'*K_test(:,i);
    end 
    
    
    %%---step 2:Perform LDA in KPCA-transformed space---%%  
    
    %%利用等式(17)求Sb%%
    y_m_all = zeros(m,1);
    y_m = zeros(m,class_num);
    kk = 0;
    for i = 1:class_num
        for j = 1:sample_num
            if train_mark(j) == 1
                kk = kk + 1;
                y_m_all = y_m_all + y_train(:,kk);
                y_m(:,i) = y_m(:,i) + y_train(:,kk);
            end
        end
        y_m(:,i) = y_m(:,i)/train_num;
    end
    y_m_all = y_m_all / train_total;
	
    Sb1 = zeros(m,class_num);
    Sb = zeros(m,m);
    for i = 1:class_num
        Sb1(:,i) = y_m(:,i) - y_m_all;
        Sb = Sb + train_num*Sb1(:,i)*Sb1(:,i)';
    end
    Sb=Sb/train_total;
    
    Sw1 = zeros(m,train_total);
    Sw = zeros(m,m);
    kk = 0;
    for i = 1:class_num
        for j = 1:sample_num
            if train_mark(j) == 1
                kk = kk + 1;
                Sw1(:,kk) = y_train(:,kk) - y_m(:,i);
                Sw = Sw + Sw1(:,kk)*Sw1(:,kk)';
            end
        end
    end
    Sw=Sw/train_total;
    
    St1 = zeros(m,train_total);
    St = zeros(m,m);
    kk = 0;
    for i = 1:class_num
        for j = 1:sample_num
            if train_mark(j) == 1
                kk = kk + 1;
                St1(:,kk) = y_train(:,kk) - y_m_all;
                St = St + St1(:,kk)*St1(:,kk)';
            end
        end
    end
    St=St/train_total;

    %clear y_m_all y_m;
    
    %%计算St的前d个最大本征值及其对应的本征向量%%
    ss = inv(Sw)*Sb;
%    ss = St;

    [V_ss, D_ss] = eig(ss);
    dd = abs(eig(ss));
    [dd_value, dd_site] = sort(dd);
    temp2 = class_num - 1; %rank(ss)
    temp1 = m - temp2 + 1;
    
    G = [];
    for i = m:-1:temp1
        v_temp = V_ss(:,dd_site(i))';
        G = [G; v_temp];
    end
    dd_bak = dd_value(temp1:m)';
    toc;
     
     %% -------- 4.calculate transform of test and train -------- %%
    orl_train_new = G * y_train;
    orl_test_new = G *y_test;

	error = [];
	for i = 1:class_num
        for j = 1:test_num
            %i,j
            kk = (i-1)*test_num + j;
            res_y = orl_test_new(:,kk);
        
            dis = [];
            for ii = 1:class_num
                for jj = 1:train_num
                    res_x = orl_train_new(:,(ii-1)*train_num+jj);
                    res = (res_y - res_x).^2;
                    temp = sqrt(sum(res));
                    res = abs(res_y - res_x);
                    temp = sum(res);
                    dis = [dis; temp];
                end
            end
            [value, site] = min(dis);
            site_new = floor((site-1)/train_num) + 1;
            
            if site_new ~= i
                error = [error; i*test_num+j];
            end
        end
	end
    right_num = [right_num; (1.0-length(error)/(class_num*test_num))*100 ];
    right_num'
end