Implement Multivariate Linear Regression with Gradient Descent and 3D Visualization in Python

Problem Background

We have a dataset of delivery times where the columns represent independent variables such as mileage and number of deliveries, and the target variable is the delivery time we want to predict. Implementing a multivariate linear regression model allows us to train on this data and make predictions.

Code Implementation

Import Required Libraries

import numpy as np
from numpy import genfromtxt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Load and Preprocess Data

data = genfromtxt(r"Delivery.csv", delimiter=",")
print(data)

x_data = data[:, :-1]  # independent variables (miles, deliveries)
print(x_data)

y_data = data[:, -1]   # target variable (delivery time)
print(y_data)

The genfromtxt function loads the CSV file into a 2‑D array where the first two columns are features and the last column is the target. x_data contains all feature columns. y_data contains the target values.

Define Model Parameters

lr = 0.0001      # learning rate
theta0 = 0       # bias term initial value
theta1 = 0       # coefficient for first feature
theta2 = 0       # coefficient for second feature
epochs = 1000    # number of gradient‑descent iterations

Learning rate (lr) controls the step size of parameter updates.

theta0, theta1, theta2 are the parameters to be optimized.

epochs determines how many times the gradient descent loop runs.

Compute Error Function

def compute_error(theta0, theta1, theta2, x_data, y_data):
    totalError = 0
    for i in range(0, len(x_data)):
        # squared error between prediction and actual value
        totalError += (y_data[i] - (theta1 * x_data[i, 0] + theta2 * x_data[i, 1] + theta0)) ** 2
    return totalError / float(len(x_data))

This function returns the mean‑squared error (MSE) of the model, which measures prediction accuracy.

Gradient Descent Algorithm

def gradient_descent_runner(x_data, y_data, theta0, theta1, theta2, lr, epochs):
    m = float(len(x_data))  # number of samples
    for i in range(epochs):
        theta0_grad = 0
        theta1_grad = 0
        theta2_grad = 0
        for j in range(0, len(x_data)):
            # compute gradients for each parameter
            theta0_grad += (1/m) * (theta1 * x_data[j, 0] + theta2 * x_data[j, 1] + theta0 - y_data[j])
            theta1_grad += (1/m) * x_data[j, 0] * (theta1 * x_data[j, 0] + theta2 * x_data[j, 1] + theta0 - y_data[j])
            theta2_grad += (1/m) * x_data[j, 1] * (theta1 * x_data[j, 0] + theta2 * x_data[j, 1] + theta0 - y_data[j])
        # update parameters using the learning rate
        theta0 = theta0 - (lr * theta0_grad)
        theta1 = theta1 - (lr * theta1_grad)
        theta2 = theta2 - (lr * theta2_grad)
    return theta0, theta1, theta2

Gradient descent iteratively adjusts theta0, theta1, and theta2 to minimize the error function until convergence.

Start Training

print("Starting theta0 = {0}, theta1 = {1}, theta2 = {2}, error = {3}".format(theta0, theta1, theta2, compute_error(theta0, theta1, theta2, x_data, y_data)))
print("Running...")

theta0, theta1, theta2 = gradient_descent_runner(x_data, y_data, theta0, theta1, theta2, lr, epochs)

print("After {0} iterations theta0 = {1}, theta1 = {2}, theta2 = {3}, error = {4}".format(epochs, theta0, theta1, theta2, compute_error(theta0, theta1, theta2, x_data, y_data)))

The script prints parameter values and error before and after training.

Initial error is high; as gradient descent runs, error decreases and parameters converge to values that better fit the data.

Visualize Model Results

ax = plt.figure().add_subplot(111, projection='3d')
ax.scatter(x_data[:, 0], x_data[:, 1], y_data, c='r', marker='o', s=100)  # data points
x0 = x_data[:, 0]
x1 = x_data[:, 1]

x0, x1 = np.meshgrid(x0, x1)
z = theta0 + x0 * theta1 + x1 * theta2  # regression plane

ax.plot_surface(x0, x1, z)  # draw regression plane
ax.set_xlabel('Miles')
ax.set_ylabel('Num of Deliveries')
ax.set_zlabel('Time')
plt.show()

A 3‑D scatter plot shows the actual data points in red.

The blue surface represents the fitted linear regression plane computed from the optimized parameters.

Conclusion

By implementing multivariate linear regression and optimizing its parameters with gradient descent, we can effectively predict delivery time from mileage and delivery count. Model accuracy depends on the dataset quality and the chosen learning rate. The 3‑D visualization provides an intuitive view of how well the model fits the data.

This end‑to‑end example—from data loading to training and visualization—demonstrates a reusable workflow for similar regression problems.