ordinary_least_squares_regressor.cpp File Reference

Linear regression example using Ordinary least squares More...

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

Include dependency graph for ordinary_least_squares_regressor.cpp:

Go to the source code of this file.

Functions
template<typename T >
std::ostream &	operator<< (std::ostream &out, std::vector< std::vector< T > > const &v)
template<typename T >
std::ostream &	operator<< (std::ostream &out, std::vector< T > const &v)
template<typename T >
bool	is_square (std::vector< std::vector< T > > const &A)
template<typename T >
std::vector< std::vector< T > >	operator* (std::vector< std::vector< T > > const &A, std::vector< std::vector< T > > const &B)
template<typename T >
std::vector< T >	operator* (std::vector< std::vector< T > > const &A, std::vector< T > const &B)
template<typename T >
std::vector< float >	operator* (float const scalar, std::vector< T > const &A)
template<typename T >
std::vector< float >	operator* (std::vector< T > const &A, float const scalar)
template<typename T >
std::vector< float >	operator/ (std::vector< T > const &A, float const scalar)
template<typename T >
std::vector< T >	operator- (std::vector< T > const &A, std::vector< T > const &B)
template<typename T >
std::vector< T >	operator+ (std::vector< T > const &A, std::vector< T > const &B)
template<typename T >
std::vector< std::vector< float > >	get_inverse (std::vector< std::vector< T > > const &A)
template<typename T >
std::vector< std::vector< T > >	get_transpose (std::vector< std::vector< T > > const &A)
template<typename T >
std::vector< float >	fit_OLS_regressor (std::vector< std::vector< T > > const &X, std::vector< T > const &Y)
template<typename T >
std::vector< float >	predict_OLS_regressor (std::vector< std::vector< T > > const &X, std::vector< float > const &beta)
void	ols_test ()
int	main ()

Detailed Description

Linear regression example using Ordinary least squares

Program that gets the number of data samples and number of features per sample along with output per sample. It applies OLS regression to compute the regression output for additional test data samples.

Author

Krishna Vedala

Definition in file ordinary_least_squares_regressor.cpp.

Function Documentation

fit_OLS_regressor()

template<typename T >

std::vector< float > fit_OLS_regressor	(	std::vector< std::vector< T > > const &	X,
		std::vector< T > const &	Y )

Perform Ordinary Least Squares curve fit. This operation is defined as

\[\beta = \left(X^TXX^T\right)Y\]

Parameters

X	feature matrix with rows representing sample vector of features
Y	known regression value for each sample

Returns

fitted regression model polynomial coefficients

Definition at line 321 of file ordinary_least_squares_regressor.cpp.

                                                            {
    // NxF
    std::vector<std::vector<T>> X2 = X;
    for (size_t i = 0; i < X2.size(); i++) {
        // add Y-intercept -> Nx(F+1)
        X2[i].push_back(1);
    }
    // (F+1)xN
    std::vector<std::vector<T>> Xt = get_transpose(X2);
    // (F+1)x(F+1)
    std::vector<std::vector<T>> tmp = get_inverse(Xt * X2);
    // (F+1)xN
    std::vector<std::vector<float>> out = tmp * Xt;
    // cout << endl
    //      << "Projection matrix: " << X2 * out << endl;
 
    // Fx1,1    -> (F+1)^th element is the independent constant
    return out * Y;
}

get_inverse()

template<typename T >

std::vector< std::vector< float > > get_inverse

(

std::vector< std::vector< T > > const &

A

)

Get matrix inverse using Row-trasnformations. Given matrix must be a square and non-singular.

Returns

inverse matrix

Definition at line 226 of file ordinary_least_squares_regressor.cpp.

                                      {
    // Assuming A is square matrix
    size_t N = A.size();
 
    std::vector<std::vector<float>> inverse(N);
    for (size_t row = 0; row < N; row++) {
        // preallocatae a resultant identity matrix
        inverse[row] = std::vector<float>(N);
        for (size_t col = 0; col < N; col++) {
            inverse[row][col] = (row == col) ? 1.f : 0.f;
        }
    }
 
    if (!is_square(A)) {
        std::cerr << "A must be a square matrix!" << std::endl;
        return inverse;
    }
 
    // preallocatae a temporary matrix identical to A
    std::vector<std::vector<float>> temp(N);
    for (size_t row = 0; row < N; row++) {
        std::vector<float> v(N);
        for (size_t col = 0; col < N; col++) {
            v[col] = static_cast<float>(A[row][col]);
        }
        temp[row] = v;
    }
 
    // start transformations
    for (size_t row = 0; row < N; row++) {
        for (size_t row2 = row; row2 < N && temp[row][row] == 0; row2++) {
            // this to ensure diagonal elements are not 0
            temp[row] = temp[row] + temp[row2];
            inverse[row] = inverse[row] + inverse[row2];
        }
 
        for (size_t col2 = row; col2 < N && temp[row][row] == 0; col2++) {
            // this to further ensure diagonal elements are not 0
            for (size_t row2 = 0; row2 < N; row2++) {
                temp[row2][row] = temp[row2][row] + temp[row2][col2];
                inverse[row2][row] = inverse[row2][row] + inverse[row2][col2];
            }
        }
 
        if (temp[row][row] == 0) {
            // Probably a low-rank matrix and hence singular
            std::cerr << "Low-rank matrix, no inverse!" << std::endl;
            return inverse;
        }
 
        // set diagonal to 1
        auto divisor = static_cast<float>(temp[row][row]);
        temp[row] = temp[row] / divisor;
        inverse[row] = inverse[row] / divisor;
        // Row transformations
        for (size_t row2 = 0; row2 < N; row2++) {
            if (row2 == row) {
                continue;
            }
            float factor = temp[row2][row];
            temp[row2] = temp[row2] - factor * temp[row];
            inverse[row2] = inverse[row2] - factor * inverse[row];
        }
    }
 
    return inverse;
}

get_transpose()

template<typename T >

std::vector< std::vector< T > > get_transpose

(

std::vector< std::vector< T > > const &

A

)

matrix transpose

Returns

resultant matrix

Definition at line 300 of file ordinary_least_squares_regressor.cpp.

                                      {
    std::vector<std::vector<T>> result(A[0].size());
 
    for (size_t row = 0; row < A[0].size(); row++) {
        std::vector<T> v(A.size());
        for (size_t col = 0; col < A.size(); col++) v[col] = A[col][row];
 
        result[row] = v;
    }
    return result;
}

is_square()

template<typename T >

bool is_square ( std::vector< std::vector< T > > const & A )

inline

function to check if given matrix is a square matrix

Returns

1 if true, 0 if false

Definition at line 59 of file ordinary_least_squares_regressor.cpp.

                                                        {
    // Assuming A is square matrix
    size_t N = A.size();
    for (size_t i = 0; i < N; i++) {
        if (A[i].size() != N) {
            return false;
        }
    }
    return true;
}

main()

int main

(

void

)

main function

Definition at line 423 of file ordinary_least_squares_regressor.cpp.

           {
    ols_test();
 
    size_t N = 0, F = 0;
 
    std::cout << "Enter number of features: ";
    // number of features = columns
    std::cin >> F;
    std::cout << "Enter number of samples: ";
    // number of samples = rows
    std::cin >> N;
 
    std::vector<std::vector<float>> data(N);
    std::vector<float> Y(N);
 
    std::cout
        << "Enter training data. Per sample, provide features and one output."
        << std::endl;
 
    for (size_t rows = 0; rows < N; rows++) {
        std::vector<float> v(F);
        std::cout << "Sample# " << rows + 1 << ": ";
        for (size_t cols = 0; cols < F; cols++) {
            // get the F features
            std::cin >> v[cols];
        }
        data[rows] = v;
        // get the corresponding output
        std::cin >> Y[rows];
    }
 
    std::vector<float> beta = fit_OLS_regressor(data, Y);
    std::cout << std::endl << std::endl << "beta:" << beta << std::endl;
 
    size_t T = 0;
    std::cout << "Enter number of test samples: ";
    // number of test sample inputs
    std::cin >> T;
    std::vector<std::vector<float>> data2(T);
    // vector<float> Y2(T);
 
    for (size_t rows = 0; rows < T; rows++) {
        std::cout << "Sample# " << rows + 1 << ": ";
        std::vector<float> v(F);
        for (size_t cols = 0; cols < F; cols++) std::cin >> v[cols];
        data2[rows] = v;
    }
 
    std::vector<float> out = predict_OLS_regressor(data2, beta);
    for (size_t rows = 0; rows < T; rows++) std::cout << out[rows] << std::endl;
 
    return 0;
}

ols_test()

void ols_test

(

)

Self test checks

Definition at line 369 of file ordinary_least_squares_regressor.cpp.

                {
    int F = 3, N = 5;
 
    /* test function = x^2 -5 */
    std::cout << "Test 1 (quadratic function)....";
    // create training data set with features = x, x^2, x^3
    std::vector<std::vector<float>> data1(
        {{-5, 25, -125}, {-1, 1, -1}, {0, 0, 0}, {1, 1, 1}, {6, 36, 216}});
    // create corresponding outputs
    std::vector<float> Y1({20, -4, -5, -4, 31});
    // perform regression modelling
    std::vector<float> beta1 = fit_OLS_regressor(data1, Y1);
    // create test data set with same features = x, x^2, x^3
    std::vector<std::vector<float>> test_data1(
        {{-2, 4, -8}, {2, 4, 8}, {-10, 100, -1000}, {10, 100, 1000}});
    // expected regression outputs
    std::vector<float> expected1({-1, -1, 95, 95});
    // predicted regression outputs
    std::vector<float> out1 = predict_OLS_regressor(test_data1, beta1);
    // compare predicted results are within +-0.01 limit of expected
    for (size_t rows = 0; rows < out1.size(); rows++) {
        assert(std::abs(out1[rows] - expected1[rows]) < 0.01);
    }
    std::cout << "passed\n";
 
    /* test function = x^3 + x^2 - 100 */
    std::cout << "Test 2 (cubic function)....";
    // create training data set with features = x, x^2, x^3
    std::vector<std::vector<float>> data2(
        {{-5, 25, -125}, {-1, 1, -1}, {0, 0, 0}, {1, 1, 1}, {6, 36, 216}});
    // create corresponding outputs
    std::vector<float> Y2({-200, -100, -100, 98, 152});
    // perform regression modelling
    std::vector<float> beta2 = fit_OLS_regressor(data2, Y2);
    // create test data set with same features = x, x^2, x^3
    std::vector<std::vector<float>> test_data2(
        {{-2, 4, -8}, {2, 4, 8}, {-10, 100, -1000}, {10, 100, 1000}});
    // expected regression outputs
    std::vector<float> expected2({-104, -88, -1000, 1000});
    // predicted regression outputs
    std::vector<float> out2 = predict_OLS_regressor(test_data2, beta2);
    // compare predicted results are within +-0.01 limit of expected
    for (size_t rows = 0; rows < out2.size(); rows++) {
        assert(std::abs(out2[rows] - expected2[rows]) < 0.01);
    }
    std::cout << "passed\n";
 
    std::cout << std::endl;  // ensure test results are displayed on screen
                             // (flush stdout)
}

operator*() [1/4]

template<typename T >

std::vector< float > operator*	(	float const	scalar,
		std::vector< T > const &	A )

pre-multiplication of a vector by a scalar

Returns

resultant vector

Definition at line 138 of file ordinary_least_squares_regressor.cpp.

                                                                      {
    // Number of rows in A
    size_t N_A = A.size();
 
    std::vector<float> result(N_A);
 
    for (size_t row = 0; row < N_A; row++) {
        result[row] += A[row] * static_cast<float>(scalar);
    }
 
    return result;
}

operator*() [2/4]

template<typename T >

std::vector< std::vector< T > > operator*	(	std::vector< std::vector< T > > const &	A,
		std::vector< std::vector< T > > const &	B )

Matrix multiplication such that if A is size (mxn) and B is of size (pxq) then the multiplication is defined only when n = p and the resultant matrix is of size (mxq)

Returns

resultant matrix

Definition at line 78 of file ordinary_least_squares_regressor.cpp.

                                                                        {
    // Number of rows in A
    size_t N_A = A.size();
    // Number of columns in B
    size_t N_B = B[0].size();
 
    std::vector<std::vector<T>> result(N_A);
 
    if (A[0].size() != B.size()) {
        std::cerr << "Number of columns in A != Number of rows in B ("
                  << A[0].size() << ", " << B.size() << ")" << std::endl;
        return result;
    }
 
    for (size_t row = 0; row < N_A; row++) {
        std::vector<T> v(N_B);
        for (size_t col = 0; col < N_B; col++) {
            v[col] = static_cast<T>(0);
            for (size_t j = 0; j < B.size(); j++) {
                v[col] += A[row][j] * B[j][col];
            }
        }
        result[row] = v;
    }
 
    return result;
}

operator*() [3/4]

template<typename T >

std::vector< T > operator*	(	std::vector< std::vector< T > > const &	A,
		std::vector< T > const &	B )

multiplication of a matrix with a column vector

Returns

resultant vector

Definition at line 112 of file ordinary_least_squares_regressor.cpp.

                                                {
    // Number of rows in A
    size_t N_A = A.size();
 
    std::vector<T> result(N_A);
 
    if (A[0].size() != B.size()) {
        std::cerr << "Number of columns in A != Number of rows in B ("
                  << A[0].size() << ", " << B.size() << ")" << std::endl;
        return result;
    }
 
    for (size_t row = 0; row < N_A; row++) {
        result[row] = static_cast<T>(0);
        for (size_t j = 0; j < B.size(); j++) result[row] += A[row][j] * B[j];
    }
 
    return result;
}

operator*() [4/4]

template<typename T >

std::vector< float > operator*	(	std::vector< T > const &	A,
		float const	scalar )

post-multiplication of a vector by a scalar

Returns

resultant vector

Definition at line 156 of file ordinary_least_squares_regressor.cpp.

                                                                      {
    // Number of rows in A
    size_t N_A = A.size();
 
    std::vector<float> result(N_A);
 
    for (size_t row = 0; row < N_A; row++) {
        result[row] = A[row] * static_cast<float>(scalar);
    }
 
    return result;
}

operator+()

template<typename T >

std::vector< T > operator+	(	std::vector< T > const &	A,
		std::vector< T > const &	B )

addition of two vectors of identical lengths

Returns

resultant vector

Definition at line 204 of file ordinary_least_squares_regressor.cpp.

                                                                     {
    // Number of rows in A
    size_t N = A.size();
 
    std::vector<T> result(N);
 
    if (B.size() != N) {
        std::cerr << "Vector dimensions shouldbe identical!" << std::endl;
        return A;
    }
 
    for (size_t row = 0; row < N; row++) result[row] = A[row] + B[row];
 
    return result;
}

operator-()

template<typename T >

std::vector< T > operator-	(	std::vector< T > const &	A,
		std::vector< T > const &	B )

subtraction of two vectors of identical lengths

Returns

resultant vector

Definition at line 183 of file ordinary_least_squares_regressor.cpp.

                                                                     {
    // Number of rows in A
    size_t N = A.size();
 
    std::vector<T> result(N);
 
    if (B.size() != N) {
        std::cerr << "Vector dimensions shouldbe identical!" << std::endl;
        return A;
    }
 
    for (size_t row = 0; row < N; row++) result[row] = A[row] - B[row];
 
    return result;
}

operator/()

template<typename T >

std::vector< float > operator/	(	std::vector< T > const &	A,
		float const	scalar )

division of a vector by a scalar

Returns

resultant vector

Definition at line 174 of file ordinary_least_squares_regressor.cpp.

                                                                      {
    return (1.f / scalar) * A;
}

operator<<() [1/2]

template<typename T >

std::ostream & operator<<	(	std::ostream &	out,
		std::vector< std::vector< T > > const &	v )

operator to print a matrix

Definition at line 22 of file ordinary_least_squares_regressor.cpp.

                                                           {
    const int width = 10;
    const char separator = ' ';
 
    for (size_t row = 0; row < v.size(); row++) {
        for (size_t col = 0; col < v[row].size(); col++) {
            out << std::left << std::setw(width) << std::setfill(separator)
                << v[row][col];
        }
        out << std::endl;
    }
 
    return out;
}

operator<<() [2/2]

template<typename T >

std::ostream & operator<<	(	std::ostream &	out,
		std::vector< T > const &	v )

operator to print a vector

Definition at line 42 of file ordinary_least_squares_regressor.cpp.

                                                               {
    const int width = 15;
    const char separator = ' ';
 
    for (size_t row = 0; row < v.size(); row++) {
        out << std::left << std::setw(width) << std::setfill(separator)
            << v[row];
    }
 
    return out;
}

predict_OLS_regressor()

template<typename T >

std::vector< float > predict_OLS_regressor	(	std::vector< std::vector< T > > const &	X,
		std::vector< float > const &	beta )

Given data and OLS model coeffficients, predict regression estimates. This operation is defined as

\[y_{\text{row}=i} = \sum_{j=\text{columns}}\beta_j\cdot X_{i,j}\]

Parameters

X	feature matrix with rows representing sample vector of features
beta	fitted regression model

Returns

vector with regression values for each sample

Definition at line 352 of file ordinary_least_squares_regressor.cpp.

  {
    std::vector<float> result(X.size());
 
    for (size_t rows = 0; rows < X.size(); rows++) {
        // -> start with constant term
        result[rows] = beta[X[0].size()];
        for (size_t cols = 0; cols < X[0].size(); cols++) {
            result[rows] += beta[cols] * X[rows][cols];
        }
    }
    // Nx1
    return result;
}