%SVO Support Vector Optimizer, low-level routine
%
% [V,J,C,NU] = SVO(K,NLAB,C,OPTIONS)
%
% INPUT
% K Similarity matrix
% NLAB Label list consisting of -1/+1
% C Scalar for weighting the errors (optional; default: 1)
% OPTIONS
% .PD_CHECK force positive definiteness of the kernel by adding a small constant
% to a kernel diagonal (default: 1)
% .BIAS_IN_ADMREG it may happen that bias of svc (b term) is not defined, then
% if BIAS_IN_ADMREG == 1, b will be taken from the midpoint of its admissible
% region, otherwise (BIAS_IN_ADMREG == 0) the situation will be considered
% as an optimization failure and treated accordingly (deafault: 1)
% .PF_ON_FAILURE if optimization is failed (or bias is undefined and BIAS_IN_ADMREG is 0)
% and PF_ON_FAILURE == 1, then Pseudo Fisher classifier will be computed,
% otherwise (PF_ON_FAILURE == 0) an error will be issued (default: 1)
%
% OUTPUT
% V Vector of weights for the support vectors
% J Index vector pointing to the support vectors
% C C which was actually used for optimization
% NU NU parameter of NUSVC algorithm, which gives the same classifier
%
% DESCRIPTION
% A low level routine that optimizes the set of support vectors for a 2-class
% classification problem based on the similarity matrix K computed from the
% training set. SVO is called directly from SVC. The labels NLAB should indicate
% the two classes by +1 and -1. Optimization is done by a quadratic programming.
% If available, the QLD function is used, otherwise an appropriate Matlab routine.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% SVC
% Copyright: D.M.J. Tax, D. de Ridder, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands
% $Id: svo.m,v 1.6 2010/02/08 15:29:48 duin Exp $
function [v,J,C,nu] = svo(K,y,C,Options)
if nargin < 4
Options = [];
end
DefOptions.pd_check = 1;
DefOptions.bias_in_admreg = 1;
DefOptions.pf_on_failure = 1;
Options = updstruct(DefOptions, Options,1);
if nargin < 3 | isempty(C)
prwarning(3,'The regularization parameter C is not specified, assuming 1.');
C = 1;
end
vmin = 1e-9; % Accuracy to determine when an object becomes the support object.
vmin1 = min(1,C)*vmin; % controls if an object is the support object
vmin2 = C*vmin; % controls if a support object is the boundary support object
% Set up the variables for the optimization.
n = size(K,1);
D = (y*y').*K;
f = -ones(1,n);
A = y';
b = 0;
lb = zeros(n,1);
ub = repmat(C,n,1);
p = rand(n,1);
D = (D+D')/2; % guarantee symmetry
% Make the kernel matrix K positive definite.
if Options.pd_check
i = -30;
while (pd_check (D + (10.0^i) * eye(n)) == 0)
i = i + 1;
end
if (i > -30),
prwarning(2,'K is not positive definite. The kernel is regularized by adding 10.0^(%d)*I',i);
end
i = i+2;
D = D + (10.0^(i)) * eye(n);
end
% Minimization procedure initialization:
% 'qp' minimizes: 0.5 x' K x + f' x
% subject to: Ax <= b
%
if (exist('qld') == 3)
v = qld (D, f, -A, b, lb, ub, p, 1);
elseif (exist('quadprog') == 2)
prwarning(1,'QLD not found, the Matlab routine QUADPROG is used instead.')
opt = optimset; opt.LargeScale='off'; opt.Display='off';
v = quadprog(D, f, [], [], A, b, lb, ub,[],opt);
else
prwarning(1,'QLD not found, the Matlab routine QP is used instead.')
verbos = 0;
negdef = 0;
normalize = 1;
v = qp(D, f, A, b, lb, ub, p, 1, verbos, negdef, normalize);
end
try
% check if the optimizer returned anything
if isempty(v)
error('Optimization did not converge.');
end
% Find all the support vectors.
J = find(v > vmin1);
Jp = J(y(J) == 1);
Jm = J(y(J) == -1);
% Sanity check: there are support objects from both classes
if isempty(J)
error('There are no support objects.');
elseif isempty(Jp)
error('There are no support objects from the positive class.');
elseif isempty(Jm)
error('There are no support objects from the negative class.');
end
% compute nu parameter
nu = sum(v(J),1)/(C*n);
% Find the SV on the boundary
I = find((v > vmin1) & (v < C-vmin2));
% Include class information into object weights
v = y.*v;
% There are boundary support objects we can use them to find a bias term
if ~isempty(I)
b = mean(y(I)-K(I,J)*v(J));
elseif Options.bias_in_admreg
% There are no boundary support objects
% We try to put the bias into the middle of admissible region
% non SV
J0 = (1:n)';
J0(J) = [];
J0p = J0(y(J0) == 1);
J0m = J0(y(J0) == -1);
% Jp and Jm are all margin errors
lb = max(y([J0p;Jm]) - K([J0p;Jm],J)*v(J));
ub = min(y([J0m;Jp]) - K([J0m;Jp],J)*v(J));
if lb > ub
error('The admissible region of the bias term is empty.');
end
prwarning(2,['The bias term is undefined. The midpoint of its admissible region is used.']);
b = (lb+ub)/2;
else
error('The bias term is undefined.');
end
v = [v(J); b];
catch
err.message = '##';
lasterror(err); % avoid problems with prwaitbar
if Options.pf_on_failure
prwarning(1,[lasterr ' Pseudo-Fisher is computed instead.']);
n = size(K,1);
%v = prpinv([K ones(n,1)])*y;
v = prpinv([K ones(n,1); ones(1,n) 0])*[y; 0];
J = [1:n]';
nu = nan;
else
rethrow(lasterror);
end
end
return;