regpathsvmoneclass.m

应用于matlab支持向量机svm的工具箱

function [alphamat,alpha0vec,lambdavec,event]=regpathsvmoneclass(xapp,kernel,kerneloption,verbose,options);


nuinit=0.999;
lambdaseuil=1.1;
lambdaseuil=size(xapp,1)*0.01;
epsilon=1e-8;
DOQP=1;
lambd=1e-8;
chouia=1e-3;


event=[];
lambdavec=[];
alphamat=[];
alpha0vec=[];

%--------------------------------------
% initialisation
%--------------------------------------
[nbtrain,dim]=size(xapp);
lambda=nuinit*nbtrain;
Kapp=normalizekernel(xapp,kernel,kerneloption) + epsilon*eye(nbtrain);


%
% INITIALIZATION
if DOQP
    c=zeros(nbtrain,1);
    A=ones(nbtrain,1);
    b=lambda;
    C= 1;
    
    indrand=randperm(nbtrain);
    Elambda=floor(lambda);
    alphainit(indrand(1:Elambda),1)=1;alphainit(indrand(Elambda+1),1)=lambda-Elambda;
    [alphaaux, multiplier, pos]=monqp(Kapp/lambda,c,A,b,C,lambd,verbose,[],[],alphainit);
    if sum(alphaaux>1 | alphaaux < 0 )>0
        fprintf('Error Init')
    end;
    alpha=zeros(nbtrain,1);
    alpha(pos)=alphaaux;
    alpha0=-lambda*multiplier;
    fx=(Kapp*alpha-alpha0)/lambda;
    lambdavec=[lambdavec lambda];
    alphamat=[alphamat alpha];
    alpha0vec=[alpha0vec alpha0];
    elbow=find(abs(fx)<lambd);
    left=find(fx<-lambd);
    right=find(fx>lambd);
else
    %--------------------------------------
    %% process fx for the parzen windows
    alpha=ones(nbtrain,1);
    alpha0=sqrt(alpha'*Kapp*alpha);
    fx=(Kapp*alpha-alpha0)/nbtrain; 
    
    indaux=find(fx<0);
    [aux,indaux1]=max(fx(indaux));
    indaux1=find(fx==aux); % see if there is any doublons
    if length(indaux1)>1
        warning('Doublons');
        indaux1=indaux1(1);
    end;
    alpha(indaux(indaux1))=lambda-floor(lambda);
    alpha0=Kapp(indaux(indaux1(1)),:)*alpha;
    fx=(Kapp*alpha-alpha0)/lambda;
    lambdavec=[lambdavec lambda];
    alphamat=[alphamat alpha];
    alpha0vec=[alpha0vec alpha0];
    elbow=[indaux(indaux1)];
    right=[];
    left=setdiff(1:nbtrain,elbow)';
end;


%
%           Main algorithm
%
while lambda > lambdaseuil & lambda > 1
    if verbose ==1
        fprintf('\rlambda=%3.3f',lambda);
    end;
    
    nbelbow=length(elbow);
    if nbelbow ~=0
        Una=[zeros(nbelbow,1);1];
        A=[Kapp(elbow,elbow) -ones(nbelbow,1); ones(1,nbelbow) 0];
        ba=A\Una;
        balpha=ba(1:end-1);
        bo=ba(end);
        
        % case elbow to 0 or 1
        lambda1=lambda + (ones(nbelbow,1)-alpha(elbow))./balpha;
        lambda2=lambda + (zeros(nbelbow,1)-alpha(elbow))./balpha;
        %case L and R
        %hellxi=sum(Kapp(:,elbow).*(ones(nbtrain,1)*(balpha)'),2) - bo;
        hellxi=Kapp(:,elbow)*balpha - bo;
        
        lambda3=lambda*(fx(left)-hellxi(left))./(-hellxi(left));
        lambda4=lambda*(fx(right)-hellxi(right))./(-hellxi(right));
        ind1=find( (lambda1+epsilon)<lambda);
        ind2=find((lambda2+epsilon)<lambda);
        ind3=find((lambda3+epsilon)<lambda);
        ind4=find((lambda4+epsilon)<lambda);     
        
        [lambdanew]=max([lambda1(ind1);lambda2(ind2);lambda3(ind3);lambda4(ind4)]);
        if isempty(lambdanew)
                fprintf('Exit due to no new lambda ... \n');
                break    
        end;
        fx=lambda/lambdanew*(fx-hellxi)+hellxi;
        alpha(elbow)=alpha(elbow)-(lambda-lambdanew)*balpha;
        alpha0=alpha0-(lambda-lambdanew)*bo;
        lambda=lambdanew;
        
        % mise