image thumbnail

Fast Linear binary SVM classifier

by

 

04 Nov 2011 (Updated )

Fast implementation of Linear binary SVM via BLAS/OpenMP API

test_lsvm

Contents

% Train binary Linear SVM with Pegasos algorithm
%
% MIN_w lambda/2 |w|^2 + 1/N SUM_i LOSS(w, X(:,i), y(i))
% where LOSS(w,x,y) = MAX(0, 1 - y w'x) is the hinge loss
%
%
% Usage
% ------
%
% w = pegasos_train(X , y , [options] );
%
%
% Inputs
% -------
%
% X                              Input data (d x N) in single/double format.
% y                              Binary label (1 x N) where y_i ={-1,1}, i=1,...,N in single/double format.
% options
%         lambda                 Regularizer (default lambda = 1/N).
%         B                      Bias term (default B = 1.0).
%         nbite                  Number of iteration (default nbite = 20*N).
%         reguperiod             Period regularization (default reguperiod = 10).
%         seed                   Seed number for internal random generator (default random seed according to time)
%
% If compiled with the "OMP" compilation flag
%         num_threads            Number of threads. If num_threads = -1, num_threads = number of core  (default num_threads = -1)
%
%


% Train binary Linear SVM classifier with coordinate descent algorithm on dual form
%
% MIN_w lambda/2 |w|^2 + 1/N SUM_i LOSS(w, X(:,i), y(i))
% where LOSS(w,x,y) = MAX(0, 1 - y w'x)² is the square hinge loss
%
%
% Usage
% ------
%
% w = cddcsvm_train(X , y , [options] );
%
%
% Inputs
% ------
%
% X                              Input data (d x N) in single/double format.
% y                              Binary label vector (1 x N) where y_i ={-1,1}, i=1,...,N in single/double format
% options
%         c                      Regularizer (default c = 1.0);
%         B                      Bias term (default B   = 0.0);
%         l1loss                 Equal to 1 if L1 loss is used instead of l2 loss (default l1loss = 0)
%         max_ite                Maximum number of iteration (default max_ite = 1000);
%         wp                     Weight for positives (default wp = 1)
%         wn                     Weight for negatives (default wn = 1)
%         eps                    Small value for convergence (default eps = 0.1)
%         tolPG                  Small value for convergence (default tolPG = 1.0e-12)
%         seed                   Seed number for internal random generator (default random seed according to time)
%
% If compiled with the "OMP" compilation flag
%
%         num_threads            Number of threads. If num_threads = -1, num_threads = number of core  (default num_threads = -1)

First example : Separate 2 Gaussians

clear, clc, close all hidden

% generate data

d                                     = 2;
Ntrain                                = 100;
m                                     = 2;

M0                                    = [0 ; 0];
R0                                    = [1 0 ; 0 1];
M1                                    = [2 ; 3];
R1                                    = [0.5 0.1 ; 0.2 1];
vect_test                             = (-4:0.1:8);

Xtrain                                = [M0(: , ones(1 , Ntrain/2)) + chol(R0)'*randn(d , Ntrain/2) , M1(: , ones(1 , Ntrain/2)) + chol(R1)'*randn(d , Ntrain/2)];
ytrain                                = [-ones(1 , Ntrain/2) , ones(1 , Ntrain/2)];

[X , Y]                               = meshgrid(vect_test);
Xtest                                 = [X(:)' ; Y(:)'];

indtrain0                             = (ytrain == -1);
indtrain1                             = (ytrain == 1);

% common option for both SVM

options.c                             = 1;
options.lambda                        = 1/(options.c*Ntrain);
options.B                             = 1;
options.l1loss                        = 0;
options.max_ite                       = 1000;
options.wp                            = 1;
options.wn                            = 1;
options.eps                           = 0.1;
options.tolPG                         = 1.0e-12;
options.nbite                         = 30*Ntrain;
options.reguperiod                    = 10;
options.seed                          = 1234543;
options.num_threads                   = -1;

% PEGASOS training

tic,wtrain                            = pegasos_train(Xtrain , ytrain , options);,toc


if(options.B~=0)
    btrain                            = wtrain(d+1);
    wtrain                            = wtrain(1:d);
    ytrain_est                        = sign(wtrain'*Xtrain + btrain);
    ytest_est                         = sign(wtrain'*Xtest  + btrain);
else
    ytrain_est                        = sign(wtrain'*Xtrain);
    ytest_est                         = sign(wtrain'*Xtest);
end

figure(1)
imagesc(vect_test , vect_test , reshape(ytest_est , length(vect_test) , length(vect_test)) )
axis ij
hold on
h = plot(Xtrain(1 , indtrain0) , Xtrain(2 , indtrain0) , 'k+' , Xtrain(1 , indtrain1) , Xtrain(2 , indtrain1) , 'm+' );
hold off
title(sprintf('PEGASOS with C = %2.2f',options.c) , 'fontsize' , 12);
legend(h , 'Negatives' , 'Positives');

% LIBLINEAR training

tic,wtrain                            = cddcsvm_train(Xtrain , ytrain , options);,toc

if(options.B~=0)
    btrain                            = wtrain(d+1);
    wtrain                            = wtrain(1:d);
    ytrain_est                        = sign(wtrain'*Xtrain + btrain);
    ytest_est                         = sign(wtrain'*Xtest  + btrain);
else
    ytrain_est                        = sign(wtrain'*Xtrain);
    ytest_est                         = sign(wtrain'*Xtest);
end

figure(2)
imagesc(vect_test , vect_test , reshape(ytest_est , length(vect_test) , length(vect_test)) )
axis ij
hold on
h = plot(Xtrain(1 , indtrain0) , Xtrain(2 , indtrain0) , 'k+' , Xtrain(1 , indtrain1) , Xtrain(2 , indtrain1) , 'm+' );
hold off
title(sprintf('LIBLINEAR with C = %2.2f',options.c) , 'fontsize' , 12);
legend(h , 'Negatives' , 'Positives');
Elapsed time is 0.004154 seconds.
Elapsed time is 0.001493 seconds.

Second example : ionosphere binary dataset

load ionosphere
[d , N]                            = size(X);

%Cross-Validation parameters

options.method                     = 7;
options.holding.rho                = 0.7;
options.holding.K                  = 10;
[Itrain , Itest]                   = sampling(X , y , options);
Ntrain                             = size(Itrain , 2);
Ntest                              = size(Itest , 2);

%SVM parameters

options.c                          = 1;
options.lambda                     = 1/(options.c*Ntrain);
options.B                          = 1;
options.l1loss                     = 0;
options.max_ite                    = 1000;
options.wp                         = 1;
options.wn                         = 1;
options.eps                        = 0.1;
options.tolPG                      = 1.0e-12;
options.nbite                      = 30*Ntrain;
options.reguperiod                 = 10;
options.seed                       = 1234543;
options.num_threads                = -1;

Acc1                               = zeros(1 , options.holding.K);
Acc2                               = zeros(1 , options.holding.K);
tp1                                = zeros(options.holding.K , 100);
fp1                                = zeros(options.holding.K , 100);
tp2                                = zeros(options.holding.K , 100);
fp2                                = zeros(options.holding.K , 100);

for cv = 1:options.holding.K

    [Xtrain , ytrain , Xtest , ytest]   = samplingset(X , y , Itrain , Itest , cv);
    wtrain1                             = pegasos_train(Xtrain , ytrain , options);
    wtrain2                             = cddcsvm_train(Xtrain , ytrain , options);

    if(options.B~=0)
        btrain1                            = wtrain1(d+1);
        wtrain1                            = wtrain1(1:d);
        fxtrain1                           = wtrain1'*Xtrain + btrain1;
        fxtest1                            = wtrain1'*Xtest  + btrain1;
        ytrain_est1                        = sign(fxtrain1);
        ytest_est1                         = sign(fxtest1);
        [tp1(cv , :) , fp1(cv , :)]        = basicroc(ytest , fxtest1);

        btrain2                            = wtrain2(d+1);
        wtrain2                            = wtrain2(1:d);
        fxtrain2                           = wtrain2'*Xtrain + btrain2;
        fxtest2                            = wtrain2'*Xtest  + btrain2;
        ytrain_est2                        = sign(fxtrain2);
        ytest_est2                         = sign(fxtest2);

        [tp2(cv , :) , fp2(cv , :)]        = basicroc(ytest , fxtest2);

    else
        fxtrain1                           = wtrain1'*Xtrain;
        fxtest1                            = wtrain1'*Xtest;
        ytrain_est1                        = sign(fxtrain1);
        ytest_est1                         = sign(fxtest1);

        fxtrain2                           = wtrain2'*Xtrain;
        fxtest2                            = wtrain2'*Xtest;
        ytrain_est2                        = sign(fxtrain2);
        ytest_est2                         = sign(fxtest2);

        [tp2(cv , :) , fp2(cv , :)]        = basicroc(ytest , fxtest2);
    end

    Acc1(cv)                              = sum(ytest_est1 == ytest)/Ntest;
    Acc2(cv)                              = sum(ytest_est2 == ytest)/Ntest;

end

mAcc1 = mean(Acc1);
mAcc2 = mean(Acc2);
mtp1  = mean(tp1);
mfp1  = mean(fp1);
mtp2  = mean(tp2);
mfp2  = mean(fp2);

fprintf('Acc PEGASOS = %2.4f , Acc LIBLINEAR = %2.4f\n', mAcc1 , mAcc2);

figure(3)
h = plot(mfp1 , mtp1 , 'k' , mfp2 , mtp2 , 'r' , 'linewidth' , 2);
grid on
axis([-0.05 , 1.05 , -0.05 , 1.05])
xlabel('False alarm rate' , 'fontsize' , 11);
ylabel('True detection rate' , 'fontsize' , 11);
title(sprintf('SVM with C = %2.2f',options.c) , 'fontsize' , 12);
legend(h , 'PEGASOS' , 'LIBLINEAR' , 'location' , 'southeast');
Acc PEGASOS = 0.8590 , Acc LIBLINEAR = 0.8648

Third example : Big data in single precision

d                                  = 1024;
N                                  = 10000;

X                                  = randn(d , N , 'single');
y                                  = single((rand(1 , N) >  0.5));
y(y==0)                            = -1;


options.c                          = 0.1;
options.lambda                     = 1/(options.c*N);
options.B                          = 1;
options.l1loss                     = 0;
options.max_ite                    = 1000;
options.wp                         = 1;
options.wn                         = 1;
options.eps                        = 0.1;
options.tolPG                      = 1.0e-12;
options.nbite                      = 30*N;
options.reguperiod                 = 10;
options.seed                       = 1234543;
options.num_threads                = -1;


tic,w1                             = pegasos_train(X , y , options);,toc
tic,w2                             = cddcsvm_train(X , y , options);,toc
Elapsed time is 0.933775 seconds.
Elapsed time is 25.907661 seconds.

Contact us