Linear Regression With Armadillo

linearregression.h

#ifndef LINEARREGRESSION_H
#define LINEARREGRESSION_H

#include "global.h"

class LinearRegression
{
public:
    void fit(arma::mat X, const arma::mat& y, double alpha = 0.1, ulong numIterations = 500) {
        prepareNormalization(X, mu, sigma);
        prepareFeatures(X);
        gradientDescent(theta,X,y,alpha,numIterations);
    }

    void fitNormal(arma::mat X, const arma::mat& y) {
        prepareNormalization(X, mu, sigma);
        prepareFeatures(X);
        normaleSolve(theta, X, y);
    }

    arma::mat predict(arma::mat X) {
        prepareFeatures(X);
        return X*theta;
    }

    double score(arma::mat X, const arma::mat& y) {
        arma::mat u = y - predict(X);
        arma::mat v = y.each_row() - arma::mean(y);
        return 1.0 - arma::accu(u.t()*u)/arma::accu(v.t()*v);
    }

private:
    void prepareFeatures(arma::mat& X) {
        // normalize entire X
        X = (X.each_row() - mu).each_row()/sigma;
        // add ones left column to X
        X = arma::join_horiz(arma::ones<arma::vec>(X.n_rows,1), X);
    }

    void gradientDescent(arma::mat& theta, const arma::mat& X, const arma::mat& y, double alpha, ulong numIterations) {
        theta.zeros(X.n_cols,y.n_cols);
        for(ulong it=1; it <= numIterations; ++it)
            theta = theta - (alpha/X.n_rows)*(X.t() * (X*theta - y));
    }

    void normaleSolve(arma::mat& theta, const arma::mat& X, const arma::mat& y) {
        theta = arma::pinv(X.t()*X)*X.t()*y;
    }

    void prepareNormalization(arma::mat& X, arma::mat& mu, arma::mat& sigma) {
        mu = arma::mean(X);
        sigma = arma::stddev(X);
    }

    arma::mat computeCost(const arma::mat& theta, const arma::mat& X, const arma::mat& y) {
        arma::mat E = X*theta - y;
        return (E.t()*E)/(2.0*X.n_rows);
    }

    arma::mat mu;
    arma::mat sigma;
    arma::vec theta;
};

#endif