MultipleLinearRegression.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef __MULTIPLE_LINEAR_REGRESSION_H__
#define __MULTIPLE_LINEAR_REGRESSION_H__
 
#include <Eigen/Eigen>
#include <cmath>
#include <vector>
 
namespace Utils
{
 
// WARNING: it doesn't let it go lower than the minimum value, if the prediction is lower, the min value is returned!
// set "TrueLinearRegression" to true if you want to use the real linear regression
class MultipleLinearRegression
{
  public:
    MultipleLinearRegression() = default;
 
    MultipleLinearRegression(const Eigen::VectorXd &initialWeights, double initialBias = 0.0, bool initialTrueLinearRegression = false)
        : W(initialWeights), b(initialBias), trueLinearRegression(initialTrueLinearRegression)
    {
    }
 
    void SetSamples(const std::vector<std::vector<double>> &x, const std::vector<double> &y)
    {
        assert(x.size() == y.size());
        assert(!x.empty());
 
        if (x.empty() || y.empty())
            return;
 
        const size_t n = std::min(x.size(), y.size());
        const size_t m = x[0].size();
 
        // if x is extended with 1 (for the bias term), then it's just a matrix thing, like this:
        // W = (X^t * X)^-1 * X^t * y
        // but from there the real W must be extracted and also b is going to be saved separately
 
        Eigen::MatrixXd X;
        X.resize(n, m + 1);
        X.col(0).setOnes();
 
        for (size_t i = 0; i < n; ++i)
            for (size_t j = 1; j <= m; ++j)
                X(i, j) = x[i][j - 1];
 
        Eigen::VectorXd Y;
        Y.resize(n);
        for (size_t i = 0; i < n; ++i)
            Y(i) = y[i];
 
        Eigen::MatrixXd Wa;
 
        /*
        Eigen::MatrixXd Xt = X.transpose();
        Eigen::MatrixXd XtX = Xt * X;
 
        if (abs(XtX.determinant()) < 1E-10)
            Wa = XtX.completeOrthogonalDecomposition().pseudoInverse() * Xt * Y;
        else
            Wa = XtX.inverse() * Xt * Y;
        */
 
        if (alpha > 0.)
        {
            // Ridge (L2) regularization: minimize ||X * w - Y||^2 + lambda * ||w||^2.
            // The penalty strength is derived from the data so it adapts to the
            // (unscaled) feature magnitudes: lambda = alpha * trace(X^t * X) / m,
            // using only the feature columns (the bias column is left unregularized).
            // It is solved by augmenting the system with sqrt(lambda) * I rows (and
            // zeros in Y), which keeps the well-conditioned QR solve instead of
            // forming X^t * X.
            double trace = 0.;
            for (size_t j = 1; j <= m; ++j)
                trace += X.col(j).squaredNorm();
 
            const double lambda = (m > 0) ? alpha * trace / static_cast<double>(m) : 0.;
 
            Eigen::MatrixXd Xa(n + m, m + 1);
            Xa.topRows(n) = X;
            Xa.bottomRows(m).setZero();
 
            const double sqrtLambda = std::sqrt(lambda);
            for (size_t j = 0; j < m; ++j)
                Xa(n + j, j + 1) = sqrtLambda;
 
            Eigen::VectorXd Ya(n + m);
            Ya.head(n) = Y;
            Ya.tail(m).setZero();
 
            Wa = Xa.colPivHouseholderQr().solve(Ya);
        }
        else
            Wa = X.colPivHouseholderQr().solve(Y);
 
        b = Wa(0);
        W.resize(m);
        for (size_t i = 0; i < m; ++i)
            W(i) = Wa(i + 1);
        // std::cout << "W: " << W << std::endl;
        // std::cout << "b: " << b << std::endl;
    }
 
    double Predict(const std::vector<double> &x) const
    {
        assert(x.size() == static_cast<size_t>(W.rows()));
 
        if (x.size() != static_cast<size_t>(W.size()))
            return 0.;
 
        Eigen::RowVectorXd xRow;
        xRow.resize(x.size());
        for (size_t i = 0; i < x.size(); ++i)
            xRow(i) = x[i];
 
        const auto val = xRow * W + b;
 
        if (trueLinearRegression)
            return val;
 
        const auto m = 1E-12;
        if (val < m)
            return m;
 
        return val;
    }
 
    void SetTrueLinearRegression(bool reg)
    {
        trueLinearRegression = reg;
    }
 
    // Ridge (L2) regularization factor (alpha). The effective penalty is derived
    // from the data as alpha * trace(X^t * X) / m, so it adapts to the (unscaled)
    // feature magnitudes. Set to 0 to disable. Applied before SetSamples.
    // A small value such as 1e-6 stabilizes collinear/rank-deficient fits without
    // materially biasing the weights.
    void SetRegularization(double reg)
    {
        alpha = reg;
    }
 
    const Eigen::VectorXd &GetWeights() const
    {
        return W;
    }
 
    double GetBias() const
    {
        return b;
    }
 
    bool IsTrueLinearRegression() const
    {
        return trueLinearRegression;
    }
 
    double GetRegularization() const
    {
        return alpha;
    }
 
  private:
    Eigen::VectorXd W;
    double b = 0.;
 
    bool trueLinearRegression = false;
    double alpha = 1e-6;
};
 
} // namespace Utils
 
#endif // __MULTIPLE_LINEAR_REGRESSION_H__