classifier_rbf_svm.hpp

Go to the documentation of this file.

#ifndef FSTCLASSIFIERRBFSVM_H
#define FSTCLASSIFIERRBFSVM_H
 /*
=========================================================================
Copyright and license:
  * FST4 software (with the exception of any bundled or externally linked
    third-party components) is Copyright (C) 2000-present Institute of
    Information Theory and Automation of the Czech Academy of Sciences (UTIA).
  * FST4 is free software: you can redistribute it and/or modify it under
    the terms of the GNU Affero General Public License as published by the
    Free Software Foundation, version 3. The full text of the license is in
    the file LICENSE in the distribution root and at
    https://www.gnu.org/licenses/agpl-3.0.html.
  * FST4 is distributed in the hope that it will be useful, but WITHOUT ANY
    WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public
    License for more details.
  * Closed-source use without the AGPL's copyleft obligations requires a
    commercial license from UTIA; in case of interest write to
    utia@utia.cas.cz. LICENSING.md in the distribution root explains both
    options, the contributor license agreement that every contribution
    requires, and the contacts (technical questions and pull requests:
    somol@utia.cas.cz, with "FST" in the subject line).
========================================================================== */
 
#include <memory>
#include <vector>
#include <cmath>
#include <limits>
#include <cassert>
#include <iostream>
#include <sstream>
#include "error.hpp"
#include "global.hpp"
#include "dense_train_tools.hpp"
#include "clonable.hpp"
#include "classifier.hpp"
#include "fst_rand.hpp"
 
/*============== Template parameter type naming conventions ==============
--------- Numeric types: -------------------------------------------------
DATATYPE - data sample values - usually real numbers (but may be integers
          in text processing etc.)
REALTYPE - must be real numbers - for representing intermediate results of
          calculations like mean, covariance etc.
IDXTYPE - index values for enumeration of data samples - (nonnegative) integers,
          extent depends on numbers of samples in data
DIMTYPE - index values for enumeration of features (dimensions), or classes (not
          class sizes) - (nonnegative) integers, usually lower extent than IDXTYPE,
          but be aware of expressions like _classes*_features*_features !
          in linearized representations of feature matrices for all classes
BINTYPE - feature selection marker type - represents ca. <10 different feature
          states (selected, deselected, sel./desel. temporarily 1st nested loop, 2nd...)
RETURNTYPE - criterion value: real value, but may be extended in future to support
          multiple values
--------- Class types: ---------------------------------------------------
SUBSET       - class of class type Subset
CLASSIFIER   - class implementing interface defined in abstract class Classifier
EVALUATOR    - class implementing interface defined in abstract class Sequential_Step
DISTANCE     - class implementing interface defined in abstract class Distance
DATAACCESSOR - class implementing interface defined in abstract class Data_Accessor
INTERVALCONTAINER - class of class type Data_Intervaller
CONTAINER    - STL container of class type Data_Interval
========================================================================== */
 
namespace FST {

template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
class Classifier_RBF_SVM : public Clonable_Impl<Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>,Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> > {
public:
        typedef Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> parent; 
        typedef Clonable_Impl<Classifier_RBF_SVM,parent> parent_impl; 
        typedef std::shared_ptr<DATAACCESSOR> PDataAccessor;      
        typedef std::shared_ptr<SUBSET> PSubset;                  
        typedef typename DATAACCESSOR::PPattern PPattern;         
        Classifier_RBF_SVM() : _C(1.0), _gamma(0.0), _l2loss(false), _max_iter(1000), _tol(0.1),
                _classes(0), _features(0), _subfeatures(0), _gamma_eff(0.0), _sv_count(0) {notify("Classifier_RBF_SVM constructor.");}
        virtual ~Classifier_RBF_SVM() {notify("Classifier_RBF_SVM destructor.");}
 
        void set_parameter_C(const REALTYPE c) {assert(c>0.0); _C=c;} 
        void set_parameter_gamma(const REALTYPE g) {assert(g>=0.0); _gamma=g;}  
        void set_l2_loss(const bool l2) {_l2loss=l2;}                
        void set_max_iter(const unsigned int mi) {assert(mi>0); _max_iter=mi;} 
        void set_tolerance(const REALTYPE t) {assert(t>0.0); _tol=t;}          
        REALTYPE get_parameter_C() const {return _C;} 
        REALTYPE get_parameter_gamma() const {return _gamma;} 
 
        void initialize(const PDataAccessor da); 

        bool optimize_parameters(const PDataAccessor da, const PSubset sub, const int max_points=100, const int max_throws=100,
                const double lgC_min=-5, const double lgC_max=9, const double lggamma_min=-15, const double lggamma_max=3, std::ostream& os=std::cout);
        bool mean_cv_accuracy(const PDataAccessor da, const PSubset sub, RETURNTYPE &mean); 

        virtual bool classify(DIMTYPE &cls, const PPattern &pattern) override;
        virtual bool classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern) override;
        virtual bool train(const PDataAccessor da, const PSubset sub) override;
        // test() is the accuracy loop inherited from the Classifier base
 
        virtual std::ostream& print(std::ostream& os) const override
        {
                os << "Classifier_RBF_SVM(C=" << _C << ", gamma=";
                if(_gamma>0.0) os << _gamma; else os << "auto(1/d)";
                os << ", loss=" << (_l2loss?"l2":"l1") << ", max_iter=" << _max_iter << ", tol=" << _tol << ")";
                return os;
        }
protected:
        // configuration
        REALTYPE _C;                  
        REALTYPE _gamma;               
        bool _l2loss;                  
        unsigned int _max_iter;        
        REALTYPE _tol;                 
 
        // data shape + trained state (support vectors only — compacted after training)
        DIMTYPE _classes;              
        DIMTYPE _features;             
        DIMTYPE _subfeatures;         
        std::vector<DIMTYPE> _index;  
        std::vector<REALTYPE> _mean;  
        std::vector<REALTYPE> _scale; 
        REALTYPE _gamma_eff;          
        IDXTYPE _sv_count;            
        std::vector<REALTYPE> _sv;    
        std::vector<REALTYPE> _beta;  
 
        // training scratch (copied by clones harmlessly)
        std::vector<REALTYPE> _Z;     
        std::vector<DIMTYPE> _y;      
        std::vector<REALTYPE> _alpha; 
        std::vector<REALTYPE> _s;     
        std::vector<REALTYPE> _K;     
        std::vector<REALTYPE> _krow;  
        std::vector<REALTYPE> _cs_x;  
        std::vector<REALTYPE> _cs_kv; 
        std::vector<RETURNTYPE> _cs_sc;      
        std::vector<REALTYPE> _beta_full; 

        static constexpr std::size_t KERNEL_BUDGET_ENTRIES=((std::size_t)1<<29)/sizeof(REALTYPE);

        REALTYPE kernel_value(const REALTYPE *a, const REALTYPE *b) const
        {
                REALTYPE dist2=0.0;
                for(std::size_t f=0;f<(std::size_t)_subfeatures;f++) {const REALTYPE dd=a[f]-b[f]; dist2+=dd*dd;}
                return std::exp(-_gamma_eff*dist2)+(REALTYPE)1.0; // +1 = the penalized-bias kernel offset
        }

        const REALTYPE* kernel_row(const std::size_t i, const std::size_t N)
        {
                if(!_K.empty()) return &_K[i*N];
                const REALTYPE *zi=&_Z[i*(std::size_t)_subfeatures];
                for(std::size_t k=0;k<N;k++) _krow[k]=kernel_value(zi,&_Z[k*(std::size_t)_subfeatures]);
                return _krow.data();
        }

        void solve_dual(REALTYPE *beta, const DIMTYPE positive_class);
        DIMTYPE model_count() const {return (_classes==2)?1:_classes;} 
};
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::initialize(const PDataAccessor da)
{
        notify("Classifier_RBF_SVM::initialize().");
        assert(da);
        assert(da->getNoOfFeatures()>0);
        if(da->getNoOfClasses()<2) throw fst_error("Classifier_RBF_SVM requires at least 2 classes");
        _classes=da->getNoOfClasses();
        _features=da->getNoOfFeatures();
        _index.resize(_features);
        _subfeatures=0; // no model yet
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::mean_cv_accuracy(const PDataAccessor da, const PSubset sub, RETURNTYPE &mean)
{
        RETURNTYPE val, result=0.0, cnt=0.0;
        for(bool b=da->getFirstSplit();b==true;b=da->getNextSplit()) {
                if(!train(da,sub)) return false;
                if(!this->test(val,da)) return false; // this-> because test() lives in the dependent Classifier base
                result+=val;
                cnt+=1.0;
        }
        mean=result/(RETURNTYPE)cnt;
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::optimize_parameters(const PDataAccessor da, const PSubset sub, const int max_points, const int max_throws, const double lgC_min, const double lgC_max, const double lggamma_min, const double lggamma_max, std::ostream& os)
{
        // the Classifier_LIBSVM::optimize_parameters protocol, C/gamma only (no coef0 axis)
        assert(da);
        assert(sub);
        assert(sub->get_d_raw()>0);
        assert(max_points>0);
        assert(max_throws>0);
        assert(lgC_max>lgC_min);
        assert(lggamma_max>lggamma_min);
        const double lgCrange=lgC_max-lgC_min;
        const double lggammarange=lggamma_max-lggamma_min;
        struct PSet {RETURNTYPE val; double Cpar, gamma;};
        std::vector<PSet> paramlist;
        // evaluate the initial (current) parameter set; gamma 0 = auto resolves to 1/d on the first train()
        RETURNTYPE val, maxval, firstval;
        if(!mean_cv_accuracy(da,sub,val)) return false;
        firstval=maxval=val;
        {
                // NOTE: this line reports _gamma_eff (auto resolved to 1/d by the train() just done);
                //       the optimizer sets _gamma explicitly from here on, so the per-candidate and
                //       "after" lines below print concrete, directly comparable values as well
                std::ostringstream sos;
                sos << "RBF SVM before optimization: accuracy=" << val << ", C=" << _C << ", gamma=" << _gamma_eff << std::endl;
                syncout::print(os,sos);
        }
        paramlist.push_back(PSet{val,(double)_C,(double)_gamma_eff});
        PSet tmpset{0.0,0.0,0.0}, bestset{0.0,0.0,0.0};
        for(int i=0;i<max_points;i++)
        { // find and evaluate next parameter set candidate (space-filling choice, as in the LibSVM twin)
                double maxdistmin=0.0;
                for(int j=0;j<max_throws;j++)
                {
                        // the bounds are natural logarithms: a candidate is exp(uniform(min,max))
                        tmpset.Cpar =std::exp(lgC_min+(((double)fst_rand()*(double)(lgC_max-lgC_min))/(double)FST_RAND_MAX));
                        tmpset.gamma=std::exp(lggamma_min+(((double)fst_rand()*(double)(lggamma_max-lggamma_min))/(double)FST_RAND_MAX));
                        double distmin=-1.0; // candidate (weighted) distance to already evaluated candidates
                        for(const PSet &ps : paramlist)
                        {
                                const double dC=(std::log(tmpset.Cpar)-std::log(ps.Cpar))/lgCrange;
                                const double dG=(std::log(tmpset.gamma)-std::log(ps.gamma))/lggammarange;
                                const double dist=std::sqrt(dC*dC+dG*dG);
                                // candidates more distant from all others are preferred; distances are
                                // weighed to penalize closeness to worse-performing candidates
                                assert(ps.val>0.0);
                                if(distmin<0 || dist*ps.val<distmin) distmin=dist*ps.val;
                        }
                        assert(distmin>=0.0);
                        if(distmin>maxdistmin) {maxdistmin=distmin; bestset=tmpset;}
                }
                // evaluate the chosen candidate and add it to the list
                set_parameter_C((REALTYPE)bestset.Cpar);
                set_parameter_gamma((REALTYPE)bestset.gamma);
                if(!mean_cv_accuracy(da,sub,val)) return false;
                bestset.val=val;
                paramlist.push_back(bestset);
                {
                        std::ostringstream sos;
                        sos << i+1 << ". ";
                        if(val>maxval) {maxval=val; sos << "BEST accuracy=";} else sos << "test accuracy=";
                        sos << bestset.val << ", C=" << _C << ", gamma=" << _gamma << std::endl;
                        syncout::print(os,sos);
                }
        }
        // identify the best parameter set among all tested ones
        assert(!paramlist.empty());
        bestset=paramlist.front();
        for(const PSet &ps : paramlist) if(ps.val>bestset.val) bestset=ps;
        set_parameter_C((REALTYPE)bestset.Cpar);
        set_parameter_gamma((REALTYPE)bestset.gamma);
        {
                std::ostringstream sos;
                sos << "RBF SVM after optimization: accuracy=" << bestset.val << ", C=" << _C << ", gamma=" << _gamma << std::endl;
                sos << "Accuracy difference=" << bestset.val-firstval << std::endl << std::endl;
                syncout::print(os,sos);
        }
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::solve_dual(REALTYPE *beta, const DIMTYPE positive_class)
{
        // kernel dual coordinate descent: min 1/2 a^T Qbar a - e^T a, 0<=a_i<=U,
        // Qbar_ij = y_i y_j K'_ij (+ diag 1/(2C) for L2-loss); the margin sums _s are
        // maintained so the gradient check is O(1); an alpha update touches one kernel row
        const std::size_t N=_y.size();
        const REALTYPE U=_l2loss?std::numeric_limits<REALTYPE>::infinity():_C;
        const REALTYPE Dii=_l2loss?(REALTYPE)1.0/((REALTYPE)2.0*_C):(REALTYPE)0.0;
        const REALTYPE Qii=(REALTYPE)2.0+Dii; // K'(x,x) = exp(0)+1 = 2 on the diagonal, always
 
        _alpha.assign(N,0.0);
        _s.assign(N,0.0);
        for(unsigned int sweep=0;sweep<_max_iter;sweep++)
        {
                REALTYPE pg_max=0.0;
                for(std::size_t i=0;i<N;i++) // FIXED cyclic order - deterministic by design (no RNG)
                {
                        const REALTYPE yi=(_y[i]==positive_class)?(REALTYPE)1.0:(REALTYPE)-1.0;
                        const REALTYPE g=yi*_s[i]-(REALTYPE)1.0+Dii*_alpha[i];
 
                        REALTYPE pg=g; // projected gradient at the box bounds
                        if(_alpha[i]<=0.0 && g>0.0) pg=0.0;
                        else if(_alpha[i]>=U && g<0.0) pg=0.0;
                        const REALTYPE apg=std::fabs(pg);
                        if(apg>pg_max) pg_max=apg;
                        if(apg>0.0)
                        {
                                const REALTYPE aold=_alpha[i];
                                REALTYPE anew=aold-g/Qii;
                                if(anew<0.0) anew=0.0; else if(anew>U) anew=U;
                                const REALTYPE step=(anew-aold)*yi;
                                _alpha[i]=anew;
                                if(step!=0.0)
                                {
                                        const REALTYPE *row=kernel_row(i,N);
                                        for(std::size_t k=0;k<N;k++) _s[k]+=step*row[k];
                                }
                        }
                }
                if(pg_max<_tol) break;
        }
        for(std::size_t i=0;i<N;i++) beta[i]=_alpha[i]*((_y[i]==positive_class)?(REALTYPE)1.0:(REALTYPE)-1.0);
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::train(const PDataAccessor da, const PSubset sub)
{
        notify("Classifier_RBF_SVM::train().");
        assert(da);
        assert(sub);
        assert(_classes==da->getNoOfClasses());
        assert(_features==da->getNoOfFeatures());
        assert(sub->get_n_raw()==da->getNoOfFeatures());
 
        // buffer the feature subset index
        _subfeatures=0;
        for(DIMTYPE fsel : sub->selected_features()) {assert(_subfeatures<_features); _index[_subfeatures++]=fsel;}
        assert(_subfeatures==sub->get_d_raw());
        assert(_subfeatures>0);
        const std::size_t d=_subfeatures;
 
        // pass 1 + standardize via the shared helper (dense_train_tools.hpp); no bias
        // coordinate here (the penalized-bias kernel offset supplies it), label rides along
        _y.clear();
        const std::size_t N=standardize_design<REALTYPE,IDXTYPE>(da,_classes,d,d,_index,_Z,_mean,_scale,
                [this](const auto *pat, const DIMTYPE c){(void)pat; _y.push_back(c);});
        if(N==0) return false;
 
        _gamma_eff=(_gamma>0.0)?_gamma:(REALTYPE)1.0/(REALTYPE)d;
 
        // kernel matrix, shared by all one-vs-rest models; above budget rows are recomputed on demand
        if(N*N<=KERNEL_BUDGET_ENTRIES)
        {
                _K.resize(N*N);
                for(std::size_t a=0;a<N;a++)
                {
                        _K[a*N+a]=(REALTYPE)2.0;
                        for(std::size_t b=a+1;b<N;b++) _K[a*N+b]=_K[b*N+a]=kernel_value(&_Z[a*d],&_Z[b*d]);
                }
                _krow.clear();
        }
        else {_K.clear(); _krow.resize(N);}
 
        // one-vs-rest models (two-class data: the single model's +1 side is class index 0)
        const DIMTYPE M=model_count();
        _beta_full.assign((std::size_t)M*N,0.0);
        for(DIMTYPE m=0;m<M;m++) solve_dual(&_beta_full[(std::size_t)m*N],m);
 
        // compact to support vectors: keep samples with a nonzero coefficient in ANY model
        _sv_count=0;
        _sv.clear(); _beta.clear();
        std::vector<std::size_t> keep;
        for(std::size_t n=0;n<N;n++)
        {
                bool used=false;
                for(DIMTYPE m=0;m<M && !used;m++) if(_beta_full[(std::size_t)m*N+n]!=0.0) used=true;
                if(used) keep.push_back(n);
        }
        _sv_count=(IDXTYPE)keep.size();
        _sv.resize(keep.size()*d);
        _beta.resize((std::size_t)M*keep.size());
        for(std::size_t k=0;k<keep.size();k++)
        {
                for(std::size_t f=0;f<d;f++) _sv[k*d+f]=_Z[keep[k]*d+f];
                for(DIMTYPE m=0;m<M;m++) _beta[(std::size_t)m*keep.size()+k]=_beta_full[(std::size_t)m*N+keep[k]];
        }
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern)
{
        assert(_subfeatures>0); // train() must have succeeded
        const std::size_t d=_subfeatures, nsv=_sv_count;
        scores.resize(_classes);
        // standardize the pattern, then evaluate the kernel against each support vector ONCE (shared by all models)
        _cs_x.resize(d); _cs_kv.resize(nsv);
        for(std::size_t f=0;f<d;f++) _cs_x[f]=((REALTYPE)pattern[_index[f]]-_mean[f])/_scale[f];
        for(std::size_t k=0;k<nsv;k++) _cs_kv[k]=kernel_value(&_sv[k*d],_cs_x.data());
        const DIMTYPE M=model_count();
        for(DIMTYPE m=0;m<M;m++)
        {
                const REALTYPE *b=&_beta[(std::size_t)m*nsv];
                REALTYPE dec=0.0; // the intercept lives inside the +1 kernel offset
                for(std::size_t k=0;k<nsv;k++) dec+=b[k]*_cs_kv[k];
                scores[m]=(RETURNTYPE)dec;
        }
        if(_classes==2) scores[1]=-scores[0]; // the single binary model scores both sides
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_RBF_SVM<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify(DIMTYPE &cls, const PPattern &pattern)
{
        if(!classify_scores(_cs_sc,pattern)) return false;
        cls=0;
        for(DIMTYPE c=1;c<_classes;c++) if(_cs_sc[c]>_cs_sc[cls]) cls=c; // first maximum wins
        return true;
}
 
 
} // namespace
#endif // FSTCLASSIFIERRBFSVM_H ///:~