svmkernelold.m

surpport vector machine,matlab

function K=svmkernel(x,kernel,kerneloption,xsup,framematrix,vector,dual,quadpen);

% Usage  K=svkernel(x,kernel,kerneloption,xsup,frame,vector,dual,quadpen);
%
% Returns the scalar product of the vectors x by using the
% mapping defined by the kernel function or x and xsup
% if the matrix xsup is defined
%
% Input
% 
% x		:input vectors
% kernel 	: kernel function
%		Type								Function					Option
%		Polynomial						'poly'					Degree (<x,xsup>+1)^d
%		Homogeneous polynomial		'polyhomog'				Degree <x,xsup>^d	
%		Gaussian							'gaussian'				Bandwidth
%		Heavy Tailed RBF				'htrbf'					[a,b]   %see Chappelle 1999	
%		Mexican 1D Wavelet 			'wavelet'
%		Frame kernel					'frame'					'sin','numerical'...	
%
%  kerneloption	: scalar or vector containing the option for the kernel
% 'gaussian' : scalar gamma is identical for all coordinates
%              otherwise is a vector of length equal to the number of 
%              coordinate
% 
%					For quadratic penalisation, add C quadratic to the end of vector
%					and C linear should be infinite
%
% 'poly' : kerneloption is a scalar given the degree of the polynomial
%          or is a vector which first element is the degree of the polynomial
%           and other elements gives the bandwidth of each dimension.
%          thus the vector is of size n+1 where n is the dimension of the problem.
%
%
% xsup		: support vector
%
% ----- 1D Frame Kernel -------------------------- 
%
%  framematrix  frame elements for frame kernel
%  vector       sampling position of frame elements
%	dual 		  dual frame
%  frame,vector and dual are respectively the matrices and the vector where the frame 
%  elements have been processed. these parameters are used only in case
%
% ------------------------------------------------  
%	Quadratic Penalization of slack variables for gaussian kernel
%  quadpen=1 default 0
%
%	see also svmreg,svmclass,svmval, kernelwavelet,kernelframe
%

% O4/O6/2000 A. Rakotomamonjy

if nargin < 8
    quadpen=0;
end;
if nargin < 6
    vector=[];
    dual=[];
end;
if nargin <5
    frame=[];
end;

if nargin<4
    xsup=x;
end;
if nargin<3
    kerneloption=1;
end;
if nargin<2
    kernel='gaussian';
end;
if isempty(xsup)
    xsup=x;
end;
[n1 n2]=size(x);
[n n3]=size(xsup);
ps  =  zeros(n1,n);			% produit scalaire
switch lower(kernel)
case 'poly'
    %fprintf('Polynomial Kernel of degree %d \n',kerneloption);
    [nk,nk2]=size(kerneloption);   
    if nk>nk2
        kerneloption=kerneloption';
        nk2=nk;
    end;
    if nk2==1
        degree=kerneloption;
        var=ones(1,n2);
        quadpen=0;
    elseif nk2 ==2
        degree=kerneloption(1);
        var=ones(1,n2)*kerneloption(2);
        quadpen=0;
    elseif nk2== n2+1
        degree=kerneloption(1);
        var=kerneloption(2:n2+1);
        quadpen=0;
    elseif nk2 ==n2+2
        degree=kerneloption(1);
        var=kerneloption(2:n2+1);
        cquadratic=kerneloption(n2+2);
        
        
        
        
        matquadpenalization=1/cquadratic*eye(n1,n1);
        quadpen=1;
    end;
    
    aux=repmat(var,n,1);
    
    ps= x *(xsup.*aux.^2)';
    K =(ps+1).^degree;
    if quadpen==1
  
        K=K+matquadpenalization;
    end;
    
case 'polyhomog'
    %fprintf('Polynomial Kernel of degree %d \n',kerneloption);
    [nk,nk2]=size(kerneloption);   
    if nk>nk2
        kerneloption=kerneloption';
        nk2=nk;
    end;
    if nk2==1
        degree=kerneloption;
        var=ones(1,n2);
    else
        if nk2 ~=n2+1
            degree=kerneloption(1);
            var=ones(1,n2)*kerneloption(2);
        else
            degree=kerneloption(1);
            var=kerneloption(2:nk2);
        end;
    end;
    aux=repmat(var,n,1);
    
    ps= x *(xsup.*aux.^2)';
    K =(ps).^degree;
    
    
case 'gaussian'
    [nk,nk2]=size(kerneloption);
    if nk ~=nk2
        if nk>nk2
            kerneloption=kerneloption';
        end;
    else
        kerneloption=ones(1,n2)*kerneloption;
    end;
    %keyboard
    %--------------------------------------------------------------------%
    %                                                                    %   
    %           Dealind with Quadratic Penalisation Case :               %   
    %                     adding 1/C on diagonal for learning            %
    %                                                                    %   
    %--------------------------------------------------------------------%
    if length(kerneloption)~=n2 & length(kerneloption)~=n2+1 
        error('Number of kerneloption is not compatible with data...');
    else
        if  length(kerneloption)==n2+1 
            cquadratic=kerneloption(n2+1);
            %matquadpenalization=1/cquadratic*eye(n1,n1);
            quadpen=1;
        else
            quadpen=0;
        end; 
        kerneloption=kerneloption(1:n2);
    end;
    
    
    
    
    
    
    for i=1:n
        ps(:,i) = sum(( (x - ones(n1,1)*xsup(i,:))./ (ones(n1,1)*kerneloption) ).^2,2);
    end;
    max=50;
    K = exp(-ps/2);
    if quadpen==1
        for i=1:n1
          

            K(i,i)=K(i,i)+1/cquadratic;
        end;
    end;
case 'htrbf'    % heavy tailed RBF  %see Chappelle Paper%
    b=kerneloption(2);
    a=kerneloption(1);
    for i=1:n
        ps(:,i) = kerneloption*sum( abs((x.^a - ones(n1,1)*xsup(i,:).^a)).^b   ,2);
    end;
    
    
    K = exp(-ps);
    
case 'wavelet'
    K=kernelwavelet(x,kerneloption,xsup);     
case 'frame'
    K=kernelframe(x,kerneloption,xsup,framematrix,vector,dual);
case 'wavelet2d'
    K=wav2dkernelint(x,xsup,kerneloption);
case 'radialwavelet2d'
    K=radialwavkernel(x,xsup);    
case 'tensorwavkernelint'
    K=tensorwavkernelint(x,xsup,kerneloption);
    
end;