Gradient Descent Optimization in Linear Regression

Hello and welcome to another session on "Regression and Gradient Descent." In today's syllabus, we will construct and fit the gradient descent algorithm into a linear regression problem. Though linear regression does have a direct solution, gradient descent is essential for computational efficiency, especially when handling larger datasets or complex models.

Gradient descent is an iterative optimization algorithm for minimizing a function, usually a loss function, quantifying the disparity between predicted and actual results. The goal of gradient descent is to find the parameters that minimize the value of the loss function. Importantly, gradient descent navigates its way to the minimum of the function by moving iteratively toward the direction of the steepest descent. However, to leverage gradient descent, the target function must be differentiable.

Gradient descent derives its name from its working mechanism: taking descents along the gradient. It operates in several iterative steps as follows:

1. Choose random values for initial parameters.
2. Calculate the cost (the difference between actual and predicted value).
3. Compute the gradient (the steepest slope of the function around that point).
4. Update the parameters using the gradient.
5. Repeat steps 2 to 4 until we reach an acceptable error rate or exhaust the maximum iterations.

A vital component of gradient descent is the learning rate, which determines the size of the descent towards the optimum solution. It is important to note that if the learning rate is too high, we may overshoot the minimum, and if it's too low, the convergence to the minimum may take too long.

Let's implement it from scratch with a basic understanding of the gradient descent algorithm. We will need two functions: one for calculating the cost and another for calculating and applying the gradient to update our parameters. Moreover, we'll add an early stop mechanism that will halt computations after a predefined number of iterations.

The cost function for linear regression, also known as the Mean Squared Error (MSE) cost, measures how well our model's predictions match the actual data. The formula is:

Here's what each part means:

• $X$ is the matrix of input features (with each row representing a data point and each column a feature).
• $\theta$ is the vector of parameters (weights) we want to learn.
• $y$ is the vector of actual target values.
• $m$ is the number of training examples.
• $X \cdot \theta$ gives us the predicted values for all data points.
• $(X\cdot\theta - y_i)$ computes the difference (error) between the predicted and actual value for each data point.
• Squaring this error ensures all differences are positive and penalizes larger errors more heavily.
• Summing over all $m$ examples gives the total squared error.
• Dividing by $2m$ gives the average squared error and the factor of $1/2$ simplifies the derivative during gradient computation.

In summary, the cost function quantifies the average squared difference between the predicted and actual values, and our goal is to find the $\theta$ that minimizes this cost.

#include <Eigen/Dense>
using namespace Eigen;

float cost(const MatrixXf& X, const VectorXf& y, const VectorXf& theta) {
    int m = y.size();
    VectorXf predictions = X * theta;
    float loss = (1.0f / (2 * m)) * (predictions - y).squaredNorm();
    return loss;
}

Next, for the gradient descent function, we follow the gradient descent update rule:

Let's break down each component of this formula:

• $\theta$: This is the vector of parameters (weights) that we want to optimize.
• $\alpha$: This is the learning rate, a small positive value that controls how big a step we take in the direction of the negative gradient.
• $m$: The number of training examples.
• $X$: The matrix of input features, where each row is a data point and each column is a feature.
• $y$: The vector of actual target values.
• $X \cdot \theta$: The vector of predicted values for all data points, using the current parameters.
• $(X\cdot\theta - y)$: The vector of errors (residuals) between the predicted and actual values.
• $X^T \cdot (X\cdot\theta - y)$: This computes the sum of the errors for each feature, weighted by the feature values. It gives us the direction and magnitude by which each parameter should be adjusted to reduce the cost.
• $\frac{1}{m}$: This averages the gradient over all training examples, ensuring the update step is not too large.

In each iteration, we compute the gradient of the cost function with respect to $\theta$, and then update $\theta$ by moving it in the direction that reduces the cost. The learning rate $\alpha$ determines how large each update step is. This process is repeated until the parameters converge to values that (locally) minimize the cost function.

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

tuple<VectorXf, VectorXf, MatrixXf> gradient_descent(
    const MatrixXf& X, const VectorXf& y,
    VectorXf theta, float alpha, int iterations
) {
    int m = y.size();
    VectorXf cost_history(iterations);
    MatrixXf theta_history(iterations, theta.size());

    for (int i = 0; i < iterations; ++i) {
        VectorXf prediction = X * theta;
        VectorXf gradient = (1.0f / m) * (X.transpose() * (prediction - y));
        theta -= alpha * gradient;
        theta_history.row(i) = theta.transpose();
        cost_history(i) = cost(X, y, theta);
    }

    return {theta, cost_history, theta_history};
}

Let's apply our gradient descent function to a simple linear regression problem. The form of linear regression is:

#include <iostream>
#include <Eigen/Dense>
#include <random>
using namespace std;
using namespace Eigen;

int main() {
    // Number of data points
    int n = 100;

    // Set up random number generators for noise
    default_random_engine generator;
    normal_distribution<float> noise(0.0, 1.0);

    // Generate input feature values (X_raw) uniformly between 0 and 2
    VectorXf X_raw(n);
    for (int i = 0; i < n; ++i) {
        X_raw(i) = 2.0f * ((float) rand() / RAND_MAX);
    }

    // Generate target values (y) using the true relationship y = 4 + 3x + noise
    VectorXf y(n);
    for (int i = 0; i < n; ++i) {
        y(i) = 4 + 3 * X_raw(i) + noise(generator);
    }

    // Prepare the feature matrix X with a column of ones (for the intercept) and the feature values
    MatrixXf X(n, 2);
    X.col(0) = VectorXf::Ones(n); // Intercept term
    X.col(1) = X_raw;             // Feature values

    // Initialize theta (parameters) randomly
    VectorXf theta = VectorXf::Random(2);

    // Set the learning rate and number of iterations
    float lr = 0.01;
    int n_iter = 1000;

    // Run gradient descent to optimize theta
    auto [final_theta, cost_hist, theta_hist] = gradient_descent(X, y, theta, lr, n_iter);

    // Output the final learned parameters
    cout << "Final theta:\n" << final_theta << endl;
    return 0;
}

Explanation of the steps:

1. Data Generation:

   • We generate 100 random input values (X_raw) between 0 and 2.
   • The target values (y) are created using the linear relationship $y = 4 + 3x + \text{noise}$, where noise is sampled from a normal distribution to simulate real-world data.

2. Feature Matrix Construction:

   • We construct the feature matrix X with two columns: the first column is all ones (to account for the intercept term), and the second column contains the generated input values.

3. Parameter Initialization:

   • The parameter vector theta (which includes both the intercept and the slope) is initialized with random values.

4. Gradient Descent Setup:

   • We set the learning rate (lr) and the number of iterations (n_iter) for the gradient descent algorithm.

5. Running Gradient Descent:

   • We call the gradient_descent function, which iteratively updates theta to minimize the cost function. It returns the optimized parameters, the history of cost values, and the history of parameter values.

6. Output:

   • Finally, we print the learned parameters (final_theta), which should be close to the true values used to generate the data (intercept ≈ 4, slope ≈ 3).

This step-by-step process demonstrates how to use gradient descent to fit a linear regression model to data, starting from data generation to parameter optimization and result interpretation.

Congratulations! You have mastered implementing the gradient descent algorithm and its application to linear regression. We covered theoretical explanations, derived the math behind the cost function and the gradient descent update rule, and brought these concepts to life by coding in C++.

It is now time to practice and solidify what you have learned. In the upcoming exercises, challenge yourself with different problems and experiment with varying parameters like the learning rate. Enjoy your journey into the world of gradients!