classifier_logreg.hpp

Go to the documentation of this file.

#ifndef FSTCLASSIFIERLOGREG_H
#define FSTCLASSIFIERLOGREG_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 <cassert>
#include <cmath>
#include <iostream>
#include "error.hpp"
#include "global.hpp"
#include "classifier.hpp"
#include "dense_train_tools.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_Logistic_Regression : public Clonable_Impl<Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>,Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> > {
public:
        typedef Clonable_Impl<Classifier_Logistic_Regression,Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> > parent_impl; 
        typedef std::shared_ptr<DATAACCESSOR> PDataAccessor;      
        typedef std::shared_ptr<SUBSET> const PSubset;            
        typedef typename DATAACCESSOR::PPattern PPattern;         
        Classifier_Logistic_Regression() : _lambda(0.01), _max_iter(100), _tol(1e-6), _classes(0), _features(0), _subfeatures(0)
                {notify("Classifier_Logistic_Regression constructor.");}
        virtual ~Classifier_Logistic_Regression() {notify("Classifier_Logistic_Regression destructor.");}
 
        void set_lambda(const REALTYPE lambda) {if(lambda<0.0) throw fst_error("Classifier_Logistic_Regression: lambda must be >=0"); _lambda=lambda;} 
        void set_max_iter(const unsigned int max_iter) {if(max_iter<1) throw fst_error("Classifier_Logistic_Regression: max_iter must be >=1"); _max_iter=max_iter;} 
        void set_tolerance(const REALTYPE tol) {if(!(tol>0.0)) throw fst_error("Classifier_Logistic_Regression: tolerance must be >0"); _tol=tol;} 
        REALTYPE get_lambda() const {return _lambda;} 
        unsigned int get_max_iter() const {return _max_iter;} 
        REALTYPE get_tolerance() const {return _tol;} 

        void initialize(const PDataAccessor da);
 
        virtual bool classify(DIMTYPE &cls, const PPattern &pattern) override;  // classifies pattern, returns the respective class index
        virtual bool classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern) override; // per-class linear scores (log-odds up to a shared shift) for ranking metrics (AUC)
        virtual bool train(const PDataAccessor da, const PSubset sub) override; // learns from designated training part of data
        // test() is the accuracy loop inherited from the Classifier base
 
        DIMTYPE get_n() const {return _features;} 
        DIMTYPE get_d() const {return _subfeatures;} 
 
        virtual std::ostream& print(std::ostream& os) const override {os << "Classifier_Logistic_Regression(lambda=" << _lambda << ", max_iter=" << _max_iter << ", tol=" << _tol << ")"; return os;}
protected:
        // knobs
        REALTYPE _lambda;        
        unsigned int _max_iter;  
        REALTYPE _tol;           
 
        // data shape (set in initialize())
        DIMTYPE _classes;             
        DIMTYPE _features;            
 
        // learned model (valid after train())
        DIMTYPE _subfeatures;         
        std::vector<DIMTYPE> _index;  
        std::vector<REALTYPE> _mean;  
        std::vector<REALTYPE> _scale; 
        std::vector<REALTYPE> _W;     
 
        // scratch buffers - contents transient, capacity persists across train() calls
        std::vector<REALTYPE> _Z;     
        std::vector<DIMTYPE> _y;      
        std::vector<REALTYPE> _P;     
        std::vector<REALTYPE> _H;     
        std::vector<REALTYPE> _grad;  
        std::vector<REALTYPE> _step;  
        std::vector<REALTYPE> _Wtry;  
        std::vector<REALTYPE> _sc;    
 
private:
        REALTYPE eval_objective(const std::vector<REALTYPE> &W, const bool want_grad);
        void assemble_hessian(const std::size_t m, const REALTYPE extra_ridge);
        // Cholesky factor/solve of the Newton system live in dense_train_tools.hpp (shared with regressor_ridge)
};
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::initialize(const PDataAccessor da)
{
        notify("Classifier_Logistic_Regression::initialize().");
        assert(da);
        assert(da->getNoOfFeatures()>0);
        if(da->getNoOfClasses()<2) throw fst_error("Classifier_Logistic_Regression 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>
REALTYPE Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::eval_objective(const std::vector<REALTYPE> &W, const bool want_grad)
{
        const std::size_t d=_subfeatures, dp1=d+1, K=_classes, N=_y.size();
        assert(N>0);
        if(want_grad) {_grad.assign(K*dp1,0.0); _P.resize(N*K);}
        REALTYPE J=0.0;
        for(std::size_t i=0;i<N;i++)
        {
                const REALTYPE *z=&_Z[i*dp1];
                REALTYPE smax; // assigned on the c==0 pass below; K>=2 is guaranteed by initialize()
                for(std::size_t c=0;c<K;c++) {
                        REALTYPE s=0.0;
                        const REALTYPE *w=&W[c*dp1];
                        for(std::size_t f=0;f<dp1;f++) s+=w[f]*z[f];
                        _sc[c]=s;
                        if(c==0 || s>smax) smax=s;
                }
                REALTYPE sumexp=0.0;
                for(std::size_t c=0;c<K;c++) {_sc[c]=std::exp(_sc[c]-smax); sumexp+=_sc[c];}
                J-=std::log(_sc[_y[i]]/sumexp);
                if(want_grad)
                        for(std::size_t c=0;c<K;c++) {
                                const REALTYPE p=_sc[c]/sumexp;
                                _P[i*K+c]=p;
                                const REALTYPE r=p-(c==(std::size_t)_y[i]?1.0:0.0);
                                REALTYPE *g=&_grad[c*dp1];
                                for(std::size_t f=0;f<dp1;f++) g[f]+=r*z[f];
                        }
        }
        J/=(REALTYPE)N;
        REALTYPE penalty=0.0;
        for(std::size_t c=0;c<K;c++)
                for(std::size_t f=0;f<d;f++) penalty+=W[c*dp1+f]*W[c*dp1+f]; // biases (f==d) unpenalized
        J+=0.5*_lambda*penalty;
        if(want_grad)
                for(std::size_t c=0;c<K;c++)
                        for(std::size_t f=0;f<dp1;f++) {
                                _grad[c*dp1+f]/=(REALTYPE)N;
                                if(f<d) _grad[c*dp1+f]+=_lambda*W[c*dp1+f];
                        }
        return J;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::assemble_hessian(const std::size_t m, const REALTYPE extra_ridge)
{
        const std::size_t d=_subfeatures, dp1=d+1, K=_classes, N=_y.size();
        _H.assign(m*m,0.0);
        for(std::size_t n=0;n<N;n++)
        {
                const REALTYPE *z=&_Z[n*dp1];
                const REALTYPE *pr=&_P[n*K];
                for(std::size_t c1=0;c1<K;c1++)
                        for(std::size_t c2=0;c2<=c1;c2++)
                        {
                                const REALTYPE coef=(c1==c2)?pr[c1]*(1.0-pr[c1]):-pr[c1]*pr[c2];
                                if(coef==0.0) continue;
                                REALTYPE *hblk=&_H[(c1*dp1)*m+c2*dp1];
                                for(std::size_t f1=0;f1<dp1;f1++)
                                {
                                        const REALTYPE zc=coef*z[f1];
                                        const std::size_t fend=(c1==c2)?f1+1:dp1; // diagonal block: lower triangle only
                                        for(std::size_t f2=0;f2<fend;f2++) hblk[f1*m+f2]+=zc*z[f2];
                                }
                        }
        }
        for(std::size_t j=0;j<m*m;j++) _H[j]/=(REALTYPE)N; // only the lower triangle is meaningful, scaling all is harmless
        for(std::size_t c=0;c<K;c++)
                for(std::size_t f=0;f<dp1;f++)
                {
                        const std::size_t j=c*dp1+f;
                        // weights get the L2 term; biases a tiny fixed ridge (softmax is
                        // shift-invariant in the biases, the gradient has no component in
                        // that null direction - the ridge only keeps Cholesky defined)
                        _H[j*m+j]+=((f<d)?_lambda:(REALTYPE)1e-8)+extra_ridge;
                }
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::train(const PDataAccessor da, const PSubset sub)
{
        notify("Classifier_Logistic_Regression::train().");
        assert(da);
        assert(sub);
        assert(sub->get_frozen_mode()==false);
        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, dp1=d+1, K=_classes, m=K*dp1;
        // Hessian memory guard: m=8192 means a 512MB system - beyond it an OOM crash
        // is likelier than a useful result, so fail with advice instead
        if(m>8192) throw fst_error("Classifier_Logistic_Regression: Newton system too large (classes*(subset size+1) > 8192); use smaller subsets (e.g. a forward search bounded to a small target subset size)");
 
        // pass 1 + standardize via the shared helper (dense_train_tools.hpp); the appended
        // bias coordinate and the label ride along in the per-pattern callback
        _y.clear();
        const std::size_t N=standardize_design<REALTYPE,IDXTYPE>(da,_classes,d,dp1,_index,_Z,_mean,_scale,
                [this](const auto *pat, const DIMTYPE c){(void)pat; _Z.push_back(1.0); _y.push_back(c);});
        if(N==0) return false;
 
        // damped Newton from W=0
        _W.assign(m,0.0);
        _sc.resize(K>dp1?K:dp1);
        _step.resize(m); _Wtry.resize(m);
        REALTYPE J=eval_objective(_W,true);
        for(unsigned int iter=0;iter<_max_iter;iter++)
        {
                REALTYPE gmax=0.0;
                for(std::size_t j=0;j<m;j++) {const REALTYPE a=std::fabs(_grad[j]); if(a>gmax) gmax=a;}
                if(gmax<_tol) break;
 
                // assemble and factor the Newton system, escalating the ridge on numerical failure
                // (in-place factoring destroys _H, so each retry reassembles). A large ridge makes the
                // Hessian approach a scaled identity, i.e. the "Newton" step degrades gracefully into a
                // short gradient step - still a descent direction, which is all the line search needs.
                REALTYPE ridge=0.0;
                assemble_hessian(m,ridge);
                bool factored=cholesky_factor(_H,m);
                while(!factored && ridge<1.0)
                {
                        ridge=(ridge==0.0)?(REALTYPE)1e-6:ridge*100.0;
                        assemble_hessian(m,ridge);
                        factored=cholesky_factor(_H,m);
                }
                if(!factored) return false; // pathological even with heavy ridge
 
                cholesky_solve(_H,m,_grad.data(),true,_step.data()); // L L^T step = -grad
 
                // step halving: guarantee monotone descent. The 60 halvings are a hard floor rather than
                // a tuning knob - at t=2^-60 the trial weights differ from the current ones by less than
                // double precision can represent, so no further halving could ever improve J.
                REALTYPE t=1.0;
                bool improved=false;
                for(int half=0;half<60;half++)
                {
                        for(std::size_t j=0;j<m;j++) _Wtry[j]=_W[j]+t*_step[j];
                        const REALTYPE Jtry=eval_objective(_Wtry,false);
                        if(Jtry<J && std::isfinite(Jtry)) {improved=true; _W.swap(_Wtry); J=Jtry; break;}
                        t*=0.5;
                }
                if(!improved) break; // FP floor - treat as converged
                J=eval_objective(_W,true); // refresh gradient and posterior cache at the accepted point
        }
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify(DIMTYPE &cls, const PPattern &pattern)
{
        notify("Classifier_Logistic_Regression::classify().");
        assert(_subfeatures>0); // train() must have succeeded
        const std::size_t d=_subfeatures, dp1=d+1, K=_classes;
        assert(_W.size()==K*dp1);
        // argmax of the linear scores == argmax of softmax posteriors (no exp needed)
        REALTYPE best;
        for(std::size_t c=0;c<K;c++)
        {
                const REALTYPE *w=&_W[c*dp1];
                REALTYPE score=w[d]; // bias
                for(std::size_t f=0;f<d;f++) score+=w[f]*(((REALTYPE)pattern[_index[f]]-_mean[f])/_scale[f]);
                if(c==0 || score>best) {best=score; cls=(DIMTYPE)c;} // first maximum wins
        }
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Logistic_Regression<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern)
{
        // the classify() linear scores, exposed per class (monotone in the softmax posteriors)
        notify("Classifier_Logistic_Regression::classify_scores().");
        assert(_subfeatures>0); // train() must have succeeded
        const std::size_t d=_subfeatures, dp1=d+1, K=_classes;
        assert(_W.size()==K*dp1);
        scores.resize(K);
        for(std::size_t c=0;c<K;c++)
        {
                const REALTYPE *w=&_W[c*dp1];
                REALTYPE score=w[d]; // bias
                for(std::size_t f=0;f<d;f++) score+=w[f]*(((REALTYPE)pattern[_index[f]]-_mean[f])/_scale[f]);
                scores[c]=(RETURNTYPE)score;
        }
        return true;
}
 
 
} // namespace
#endif // FSTCLASSIFIERLOGREG_H ///:~