drwnLinearRegressor.h

Go to the documentation of this file.

/******************************************************************************
** DARWIN: A FRAMEWORK FOR MACHINE LEARNING RESEARCH AND DEVELOPMENT
** Distributed under the terms of the BSD license (see the LICENSE file)
** Copyright (c) 2007-2015, Stephen Gould
** All rights reserved.
**
******************************************************************************
** FILENAME:    drwnLinearRegressor.h
** AUTHOR(S):   Stephen Gould <stephen.gould@anu.edu.au>
**
*****************************************************************************/

#pragma once

#include <cstdlib>
#include <vector>

#include "drwnBase.h"
#include "drwnFeatureMaps.h"
#include "drwnRegression.h"
#include "drwnOptimizer.h"

using namespace std;


// drwnLinearRegressorBase --------------------------------------------------

class drwnLinearRegressorBase : public drwnRegression, protected drwnOptimizer {
 public:
    static double HUBER_BETA;    
    static double REG_STRENGTH;  
    static int MAX_ITERATIONS;   

 protected:
    VectorXd _theta;    
    int _penalty;       
    double _beta;       
    int _regularizer;   
    double _lambda;     

    // cached data for parameter estimation
    // \todo change to drwnDataset when ownership flag is implemented
    const vector<vector<double> > *_features;
    const vector<double> *_targets;
    const vector<double> *_weights;

 public:
    drwnLinearRegressorBase();
    drwnLinearRegressorBase(unsigned n);
    drwnLinearRegressorBase(const drwnLinearRegressorBase &r);
    ~drwnLinearRegressorBase();

    // access functions
    virtual const char *type() const { return "drwnLinearRegressor"; }

    // i/o
    virtual bool save(drwnXMLNode& xml) const;
    virtual bool load(drwnXMLNode& xml);

    // training
    using drwnRegression::train;
    virtual double train(const drwnRegressionDataset& dataset);

    // evaluation (regression)
    virtual double getRegression(const vector<double>& features) const = 0;

 protected:
    // drwnOptimizer interface
    double objective(const double *x) const;
    void gradient(const double *x, double *df) const;
    virtual double objectiveAndGradient(const double *x, double *df) const = 0;
};

// drwnTLinearRegressor -----------------------------------------------------

template<class FeatureMap = drwnBiasFeatureMap>
class drwnTLinearRegressor : public drwnLinearRegressorBase {
 public:
    drwnTLinearRegressor() : drwnLinearRegressorBase() { /* do nothing */ }
    drwnTLinearRegressor(unsigned n) :
        drwnLinearRegressorBase(n) { initialize(n); }
    drwnTLinearRegressor(const drwnTLinearRegressor<FeatureMap> &c) :
       drwnLinearRegressorBase(c) { /* do nothing */ }

    ~drwnTLinearRegressor() { /* do nothing */ }

    // access
    virtual drwnTLinearRegressor<FeatureMap> *clone() const {
        return new drwnTLinearRegressor<FeatureMap>(*this);
    }

    // initialization
    virtual void initialize(unsigned n);

    // evaluation (regression)
    virtual double getRegression(const vector<double>& features) const;

 protected:
    virtual double objectiveAndGradient(const double *x, double *df) const;
};

// drwnLinearRegressor ------------------------------------------------------

typedef drwnTLinearRegressor<> drwnLinearRegressor;

// drwnTLinearRegressor implementation --------------------------------------

template<class FeatureMap>
void drwnTLinearRegressor<FeatureMap>::initialize(unsigned n)
{
    drwnRegression::initialize(n);
    const FeatureMap phi(_nFeatures);
    _theta = VectorXd::Zero(phi.numParameters());
}

template<class FeatureMap>
double drwnTLinearRegressor<FeatureMap>::getRegression(const vector<double>& features) const
{
    DRWN_ASSERT(features.size() == (unsigned)_nFeatures);

    vector<double> t(_theta.rows());
    Eigen::Map<VectorXd>(&t[0], t.size()) = _theta;

    const FeatureMap phi(_nFeatures);
    return phi.dot(t, features);
}

template<class FeatureMap>
double drwnTLinearRegressor<FeatureMap>::objectiveAndGradient(const double *x, double *df) const
{
    // compute gradient and objective
    double obj = 0.0;

    const unsigned m = _targets->size();
    const FeatureMap phi(_nFeatures);
    const vector<double> vx(x, x + _n);
    vector<double> vdf(_n, 0.0);

    if (_penalty == 0) {
        // L2 penalty
        for (unsigned i = 0; i < m; i++) {
            double predicted = phi.dot(vx, (*_features)[i]);
            double dist =  predicted - (*_targets)[i];
            double wdist = (_weights == NULL) ? dist : dist * (*_weights)[i];
            obj += dist * wdist;
            phi.mac(vdf, (*_features)[i], wdist);
        }

        obj *= 0.5;
    } else {
        // huber penalty
        double dh;
        for (unsigned i = 0; i < m; i++) {
            double predicted = phi.dot(vx, (*_features)[i]);
            double u = predicted - (*_targets)[i];
            if (_weights == NULL) {
                obj += drwn::huberFunctionAndDerivative(u, &dh, _beta);
            } else {
                obj += (*_weights)[i] * drwn::huberFunctionAndDerivative(u, &dh, _beta);
                dh *= (*_weights)[i];
            }
            phi.mac(vdf, (*_features)[i], dh);
        }
    }

    memcpy((void *)df, (void *)&vdf[0], _n * sizeof(double));
    if (m == 0.0) return 0.0;

    obj /= (double)m;
    Eigen::Map<VectorXd>(df, _n) /= (double)m;

    // regularization
    switch (_regularizer) {
    case 0: // sum-of-squares
        {
            double weightNorm = 0.0;
            for (unsigned i = 0; i < _n; i++) {
                weightNorm += x[i] * x[i];
                df[i] += _lambda * x[i];
            }

            obj += 0.5 * _lambda * weightNorm;
        }
        break;

    case 1: // huber
        {
            double dh;
            for (unsigned i = 0; i < _n; i++) {
                obj += _lambda * drwn::huberFunctionAndDerivative(x[i], &dh, 1.0e-3);
                df[i] += _lambda * dh;
            }
        }
        break;

    default:
        DRWN_LOG_ERROR("unsupported regularizer " << _regularizer);
    }

    return obj;
}