Implementing Multiple Linear Regression from Scratch

Introduction

Welcome to our exciting second class in the Regression and Gradient Descent series! In the previous lesson, we covered Simple Linear Regression. Now, we're transitioning toward Multiple Linear Regression, a powerful tool for examining the relationship between a dependent variable and several independent variables.

Consider a case where we need to predict house prices, which undoubtedly depend on multiple factors, such as location, size, and the number of rooms. Multiple Linear Regression accounts for these simultaneous predictors. In today's lesson, you'll learn how to implement this concept in C++!

Multiple Linear Regression - The Concept

Multiple Linear Regression builds upon the concept of Simple Linear Regression, accounting for more than one independent variable.

Let's recall the Simple Linear Regression equation:

y=β0+β1xy = \beta_0 + \beta_1x

For Multiple Linear Regression, we add multiple independent variables, x1,x2,...xmx_1, x_2, ... x_m:

Linear Algebra Behind: Dataset Representation

Suppose we had n data points (equations), each with m features (x values) Then X would look like:

X=[1x1,1x1,2…x1,m1x2,1x2,2…x2,m⋮⋮⋮⋱⋮1xn,1xn,2…xn,m]\mathbf{X} = \begin{bmatrix} 1 & x_{1,1} & x_{1,2} & \ldots & x_{1,m} \\ 1 & x_{2,1} & x_{2,2} & \ldots & x_{2,m} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n,1} & x_{n,2} & \ldots & x_{n,m} \\ \end{bmatrix}

Each row represents the m features for a single data point. Notice how we include a column of 1's the represent the intercept (also called bias) of each equation.

For each row (equation), there is a corresponding y value. So y looks like:

y=[y1y2â‹®yn]\mathbf{y} = \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix}

The normal equation results in a vector:

[β0β1⋮βm]\begin{bmatrix} \mathbf{β}_0 \\ \mathbf{β}_1 \\ \vdots \\ \mathbf{β}_{m} \end{bmatrix}

Linear Algebra Behind: Making a Prediction

Now, for any set of features x1{x_{1}} through xm{x_{m}}, we can predict the y^\hat{y} value as:

y^=(1⋅β0)+(β1⋅x1)+(β2⋅x2)+...+(βm⋅xm)\hat{y} = (1 \cdot {β}_0) + ({β}_1 \cdot x_{1}) + ({β}_2 \cdot x_{2}) + ... + ({β}_m \cdot x_{m})

To calculate all the predictions at once, we take the dot product of X{X} and β{β}

y=[y1y2⋮yn]=[1x1,1x1,2…x1,m1x2,1x2,2…x2,m⋮⋮⋮⋱⋮1xn,1xn,2…xn,m][β0β1⋮βm]=X⋅β\mathbf{y} = \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix} = \begin{bmatrix} 1 & x_{1,1} & x_{1,2} & \ldots & x_{1,m} \\ 1 & x_{2,1} & x_{2,2} & \ldots & x_{2,m} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n,1} & x_{n,2} & \ldots & x_{n,m} \\ \end{bmatrix} \begin{bmatrix} \beta_{0} \\ \beta_{1} \\ \vdots \\ \beta_{m} \end{bmatrix} = X \cdot \mathbf{\beta}

Linear Algebra Behind: Math Solution

To implement Multiple Linear Regression, we'll leverage some Linear Algebra concepts. Using the Normal Equation, we can calculate the coefficients for our regression equation:

β=(XTX)−1XTy\beta = (X^T X)^{-1} X^T y

Where XX is a matrix of features and yy is a vector of the target variable values. Like Simple Linear Regression, residuals (the differences between actual and predicted values) play a significant role. The smaller these residuals, the better the model fits.

Implementing Multiple Linear Regression from Scratch

Let's roll up our sleeves and start coding! We'll primarily rely on Eigen to handle numerical operations and matrices.

First, we set up our dataset:

C++
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main() {
    MatrixXf X(5, 3);
    X << 73, 67, 43,
         91, 88, 64,
         87, 134, 58,
         102, 43, 37,
         69, 96, 70;

    VectorXf y(5);
    y << 56, 81, 119, 22, 103;

Next, we calculate our matrix of coefficients, β\boldsymbol{\beta}, using the Normal Equation:

  • Enhance our feature matrix, XX, with an extra column of ones to account for the intercept.
C++
    MatrixXf X_b(X.rows(), X.cols() + 1);
    X_b << VectorXf::Ones(X.rows()), X;
  • Compute the coefficients β\boldsymbol{\beta} using the Normal Equation.
C++
    VectorXf beta = (X_b.transpose() * X_b).inverse() * X_b.transpose() * y;

Model's Performance Evaluation

After completing our model, we need to evaluate its performance. We use the coefficient of determination (R2R^2 score) for this purpose. It indicates how well our model fits the data. The formula is:

R2=1−SSresidualsSStotalR^2 = 1 - \frac{SS_{\text{residuals}}}{SS_{\text{total}}}

Here, SSresidualsSS_{\text{residuals}} is the residual sum of squares, and SStotalSS_{\text{total}} is the total sum of squares:

SSresiduals=∑i=1n(yi−y^i)2SS_{\text{residuals}} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2

where yiy_i are the observed values and y^i\hat{y}_i are the predicted values from the regression model.

SStotal=∑i=1n(yi−yˉ)2SS_{\text{total}} = \sum_{i=1}^{n} (y_i - \bar{y})^2

where yiy_i are the observed values and yˉ\bar{y} is the mean of the observed data.

A higher R2R^2 value (closer to 1) indicates a better model fit.

C++
    VectorXf predictions = X_b * beta;
    float ss_residuals = (y - predictions).squaredNorm();
    float ss_total = (y.array() - y.mean()).square().sum();
    float r2_score = 1 - (ss_residuals / ss_total);

    cout << "R^2 Score: " << r2_score << endl;
    return 0;
}

Lesson Summary and Practice

Congratulations on mastering Multiple Linear Regression! You've effectively bridged the gap from concept to implementation, designing a regression model in C++ from scratch.

Prepare for the upcoming lesson to delve more deeply into Regression Analysis. Meanwhile, make sure to practice and refine your newly acquired skills!

Sign up

Join the 1M+ learners on CodeSignal

Be a part of our community of 1M+ users who develop and demonstrate their skills on CodeSignal