hesse.m

用MATLAB编写的基于拟牛顿算法的最优分类器程序

%计算目标函数的海赛矩阵
function H=hesse(A,Y,w,r,m,n,kk)
%A为sample函数产生的样本矩阵
%Y为A中每个样本的分类标号
%w和r分别为待求的n×1列向量和标量，第一次为给定的初始值，以后为迭代得到的值
%m和n分别为样本的个数及属性
%kk是来标志用哪个光滑多项式代替正号函数
D=diag(Y);
e=ones(m,1);
DA=D*A;
De=D*e;
V=e-(DA*w)+(De*r);
H=zeros(n+1,n+1);
for i=1:m
  if V(i)>=0.01
    for j=1:n
      for k=1:n
        H(j,k)=H(j,k)+D(i,i)*A(i,k)*D(i,i)*A(i,j);
        end
     end
    for s=1:n
      H(n+1,s)=H(n+1,s)-D(i,i)*A(i,s)*D(i,i);
      H(s,n+1)=H(n+1,s);
    end
    H(n+1,n+1)= H(n+1,n+1)+D(i,i)*D(i,i);
    elseif V(i)<=0.01&V(i)>=-0.01
        switch kk
            case 1 %用多项式1代替正号函数
             p1V=V(i)^2*25+V(i)*0.5+0.0025;
             p11=50*V(i)+0.5;
            for j=1:n
              for k=1:n
                H(j,k)=H(j,k)+p11*p11*D(i,i)*A(i,j)*D(i,i)*A(i,k)+50*D(i,i)*A(i,j)*D(i,i)*A(i,k);
              end
            end
            for s=1:n
             H(n+1,s)=H(n+1,s)-p11*p11*D(i,i)*D(i,i)*A(i,s)-50*p1V*D(i,i)*D(i,i)*A(i,s);
             H(s,n+1)=H(n+1,s); 
            end
            H(n+1,n+1)=H(n+1,n+1)+p11*p11*D(i,i)*D(i,i)+50*p1V*D(i,i)*D(i,i);
            case 2 %用多项式2代替正号函数
             pV2=(-1/1600)*(100*V(i)+1)^3*(100*V(i)-3);
             p21=(-3/16)*(100*V(i)+1)^2*(100*V(i)-3);
             p22=(-1/16)*(100*V(i)+1)^3;
             p23=(-600/16)*(100*V(i)+1)*(100*V(i)-3);
             p24=(-300/16)*(100*V(i)+1)^2;
             for j=1:n
               for k=1:n
               H(j,k)=H(j,k)+(p21+p22)*(p21+p22)*D(i,i)*A(i,j)*D(i,i)*A(i,k)+pV2*(p23+p24)*D(i,i)*A(i,j)*D(i,i)*A(i,k);
               end
             end
             for s=1:n
                H(n+1,s)=H(n+1,s)-(p21+p22)*(p21+p22)*D(i,i)*D(i,i)*A(i,s)-pV2*(p23+p24)*D(i,i)*D(i,i)*A(i,s);
                H(s,n+1)=H(n+1,s);
             end
             H(n+1,n+1)=H(n+1,n+1)+(p21+p22)*(p21+p22)*D(i,i)*D(i,i)+pV2*(p23+p24)*D(i,i)*D(i,i);
            case 3 %用多项式3代替正号函数
            pV3=(1/32)*(100*V(i)+1)^3*(10000*V(i)^3-300*V(i)^2+V(i)+0.05);
            p31=(300/32)*(100*V(i)+1)^2*(10000*V(i)^3-300*V(i)^2+V(i)+0.05);
            p32=(1/32)*(100*V(i)+1)^3*(30000*V(i)^2-600*V(i)+1);
            p33=(60000/32)*(100*V(i)+1)*(10000*V(i)^3-300*V(i)^2+V(i)+0.05);
            p34=(300/32)*(100*V(i)+1)^2*(30000*V(i)^2-600*V(i)+1);
            p35=(300/32)*(100*V(i)+1)^2*(30000*V(i)^2-600*V(i)+1);
            p36=(1/32)*(100*V(i)+1)^3*(60000*V(i)-600);
            for j=1:n
              for k=1:n
                H(j,k)=H(j,k)+(p31+p32)*(p31+p32)*D(i,i)*A(i,j)*D(i,i)*A(i,k)+pV3*(p33+p34+p35+p36)*D(i,i)*A(i,j)*D(i,i)*A(i,k);
              end
            end
            for s=1:n
                 H(n+1,s)=H(n+1,s)-(p31+p32)*(p31+p32)*D(i,i)*A(i,s)*D(i,i)-pV3*(p33+p34+p35+p36)*D(i,i)*A(i,s)*D(i,i);
                 H(s,n+1)=H(n+1,s);
            end
             H(n+1,n+1)= H(n+1,n+1)+(p31+p32)*(p31+p32)*D(i,i)*D(i,i)+pV3*(p33+p34+p35+p36)*D(i,i)*D(i,i);
             otherwise
                 disp('最后一个参数kk的值只能取1,2,3!');
        end
  end
end
H=H+eye(n+1);