Building Linear Regression from Scratch: A Mathematical Journey with SymPy

Understand Machine Learning Through Mathematical First Principles

Machine learning often feels like a black box. We import libraries, call fit methods, and get predictions without truly understanding what happens under the hood. In this article, we’ll demystify linear regression by building it from mathematical first principles using symbolic computation.

Why Build from Scratch?

While libraries like scikit-learn are powerful and production-ready, implementing algorithms from scratch offers invaluable insights into how they actually work. Using SymPy for symbolic mathematics allows us to see and manipulate the actual mathematical expressions, making the learning process transparent and intuitive.

The Mathematical Foundation

Linear Regression Formula

Linear regression finds the best straight line that fits through our data. The line can be described with a simple equation:

Where:

• h = The value we want to predict (salary)
• x = The input variable (years of experience)
• m = The slope (how much y changes when x increases by 1)
• c = The intercept (where the line crosses the y-axis)

The model works by finding the values of m and c that make the line fit the data as closely as possible. We do this by minimizing the total squared distance between the actual data points and the predicted line.

Cost Function (Mean Squared Error)

The cost function measures how far our predictions are from the actual values. It’s the function we want to minimize to find the best line.

Formula:

Expanded form:

How it works:

1. Calculate errors: For each data point, find the difference between the predicted value (mx + c) and the actual value (y)
2. Square the errors: Square each difference to penalize larger errors more heavily and make all errors positive
3. Sum squared errors: Add up all the squared differences
4. Average the errors: Divide by 2n to normalize the cost (the factor of 2 is for mathematical convenience in derivatives)
5. Minimize: Find the values of m and c that minimize J(m, c) to get the best-fitting line

The division by 2n ensures that the cost doesn’t grow with the size of the dataset, making it comparable across different datasets.

Implementation: Step by Step

Setting Up Our Environment

First, let’s import the necessary libraries:

import sympy as sp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

STEP 1: Define Symbolic Variables

We initialize SymPy for mathematical notation and define our symbolic variables:

sp.init_printing(use_unicode=False, use_latex="mathjax")
m, c, x, y = sp.symbols('m c x y')

STEP 2: Defining the Hypothesis in Symbolic Notation

hypothesis = m * x + c
print("\nStep 2: Hypothesis Formula")
sp.pprint(hypothesis)

This creates our linear equation symbolically, allowing us to manipulate it mathematically.

STEP 3: Defining the Cost Function (MSE)

cost_func = (hypothesis - y)**2
print("\nStep 3: Cost Function (Single Point MSE)")
sp.pprint(cost_func)

This represents the squared error for a single data point.

STEP 4: Computing Gradients with Calculus

Here’s where SymPy shines — we can compute derivatives symbolically:

# Create unevaluated derivative objects for display
grad_m_eqn = sp.Derivative(cost_func, m)
grad_c_eqn = sp.Derivative(cost_func, c)

# Calculate the actual derivative expressions
grad_m_expr = grad_m_eqn.doit()
grad_c_expr = grad_c_eqn.doit()

print("\nStep 4: Symbolic Gradients (Calculus)", end="\n")
print("Partial Derivative for Weight (m):", end="\n")

print("\n\n")
sp.pprint(grad_m_eqn)
print("Evaluates to:")
sp.pprint(grad_m_expr)

print("\nPartial Derivative for Bias (c):")

print("\n\n")
sp.pprint(grad_c_eqn)
print("Evaluates to:")
sp.pprint(grad_c_expr)

These gradients tell us how to adjust m and c to minimize our cost function.

STEP 5: Loading Data with Pandas

print("\nStep 5: Loading Data with Pandas", end="\n")
try:
    # Replace 'data.csv' with your actual file path
    df = pd.read_csv('../data/Salary_Data.csv')
    X_train = df['YearsExperience'].values
    Y_train = df['Salary'].values
    print("\nData successfully loaded from CSV using Pandas.")

except FileNotFoundError:
    print("\nCSV not found. Using placeholder data for demonstration.")
    X_train = np.array([1, 2, 3, 4, 5])
    Y_train = np.array([3.1, 4.9, 7.2, 8.8, 11.1])

print("X_train = " ,X_train)
print("Y_train = ", Y_train)

STEP 6: Gradient Descent Training

def train_linear_regression(data_x, data_y, lr=0.01, epochs=1000):
    curr_m, curr_c = 0.0, 0.0
    n = len(data_x)
    print(f"\nStep 6: Training (LR={lr}, Epochs={epochs})...")

    for epoch in range(epochs):
        sum_grad_m = 0
        sum_grad_c = 0
        for xi, yi in zip(data_x, data_y):
            subs = {m: curr_m, c: curr_c, x: xi, y: yi}
            sum_grad_m += grad_m_expr.subs(subs)
            sum_grad_c += grad_c_expr.subs(subs)

        curr_m -= lr * (sum_grad_m / n)
        curr_c -= lr * (sum_grad_c / n)

        if epoch % 250 == 0:
            print(f"  Epoch {epoch}: m = {float(curr_m):.4f}, c = {float(curr_c):.4f}")

    return float(curr_m), float(curr_c)

This function implements gradient descent:

1. We start with random values for m and c (both 0.0)
2. For each epoch, we compute the gradients for all data points
3. We update m and c in the direction that reduces the cost
4. We repeat until convergence

STEP 7: Evaluation and Visualization

def evaluate_and_plot(data_x, data_y, fm, fc):
    y_true = np.array(data_y)
    y_pred = fm * data_x + fc

    # Accuracy Metrics
    r2 = 1 - (np.sum((y_true - y_pred)**2) / np.sum((y_true - np.mean(y_true))**2))
    mae = np.mean(np.abs(y_true - y_pred))

    print("\nStep 7: Accuracy Metrics")
    print(f"  - R² Score: {r2:.4f}")
    print(f"  - Mean Absolute Error (MAE): {mae:.4f}")

    # Plotting
    plt.figure(figsize=(10, 6))
    plt.scatter(data_x, data_y, color='red', label='Actual Data')
    plt.plot(data_x, y_pred, color='blue', label=f'Model: y={fm:.2f}x + {fc:.2f}')
    plt.title('Linear Regression Fit (SymPy + Pandas)')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

STEP 8: Computing Total Cost

cost_lambda = sp.lambdify((m, c, x, y), cost_func, 'numpy')

def compute_total_cost(m_val, c_val, x_data, y_data):
    # Calculates the average MSE for a grid of m and b values
    total_error = 0
    for xi, yi in zip(x_data, y_data):
        total_error += cost_lambda(m_val, c_val, xi, yi)
    return total_error / len(x_data)

STEP 9: Contour Plot Visualization

To truly understand how gradient descent finds the optimal solution, we can visualize the cost function landscape:

def plot_cost_contour(data_x, data_y, fm, fc):
    # 1. Create a grid of m and b values around our final solution
    m_range = np.linspace(fm - 2, fm + 2, 50)
    c_range = np.linspace(fc - 2, fc + 2, 50)
    M, C = np.meshgrid(m_range, c_range)

    # 2. Compute cost for every point on the grid
    Z = np.array([compute_total_cost(mv, cv, data_x, data_y) for mv, cv in zip(np.ravel(M), np.ravel(C))])
    Z = Z.reshape(M.shape)

    # 3. Plotting
    plt.figure(figsize=(8, 6))
    cp = plt.contourf(M, C, Z, levels=20, cmap='viridis') # Filled contours
    plt.colorbar(cp, label='Cost (MSE)')

    # Mark the final optimized point
    plt.plot(fm, fc, 'ro', label=f'Minimum (m={fm:.2f}, b={fc:.2f})')

    plt.title('Cost Function Contour Landscape')
    plt.xlabel('Slope (m)')
    plt.ylabel('Intercept (c)')
    plt.legend()
    plt.show()

final_m, final_c = train_linear_regression(X_train, Y_train)
evaluate_and_plot(X_train, Y_train, final_m, final_c)
plot_cost_contour(X_train, Y_train, final_m, final_c)ape. The red dot marks where gradient descent found the minimum cost — our optimal solution.

Key Insights

1. Transparency Through Symbolic Math

Using SymPy allows us to see the actual mathematical expressions we’re working with, making the algorithm interpretable rather than a black box.

2. Understanding Gradient Descent

By implementing gradient descent from scratch, we understand how the algorithm iteratively finds the optimal parameters by following the slope of the cost function.

3. The Power of Visualization

The contour plot reveals the landscape of the cost function, showing us visually why gradient descent works — it’s literally descending down the slopes to find the valley (minimum cost).

4. From Theory to Practice

This implementation bridges theoretical machine learning concepts with practical coding skills, demonstrating that complex algorithms are built on simple mathematical principles.

Performance Considerations

While our implementation is educational and transparent, production systems typically use optimized libraries like scikit-learn because:

• They’re highly optimized in C/C++
• They handle edge cases and numerical stability
• They include advanced features like regularization
• They’re battle-tested on diverse datasets

However, understanding the fundamentals through implementations like ours makes you a better practitioner who can debug, optimize, and extend these algorithms when needed.

Applications

Linear regression, despite its simplicity, has wide applications:

• Salary Prediction (as demonstrated here)

It’s often the first step in more complex models and remains relevant even in the age of deep learning.

Conclusion

Building machine learning algorithms from scratch using symbolic mathematics provides deep insights that complement the use of high-level libraries. This approach:

• Demystifies how algorithms work
• Strengthens mathematical intuition
• Enables better debugging and optimization
• Builds confidence in applying ML to real-world problems

The complete implementation demonstrates that powerful predictive models can be built from simple mathematical principles, one gradient step at a time.

Tools & Technologies

• Python: Programming language for implementation
• SymPy: For symbolic mathematics and automatic differentiation
• NumPy: For efficient numerical computations
• Pandas: For data manipulation and analysis
• Matplotlib: For comprehensive visualization

Acknowledgments

This work was made possible through the collaboration and guidance of talented colleagues and experienced mentors. Special thanks to our mentors for their invaluable insights and support throughout this learning journey.

Have you implemented machine learning algorithms from scratch? What insights did you gain? Share your experiences in the comments!

About the Authors:

***Harish Gopalakrishnan*** is a Senior Software Engineer at Wipro with expertise in software development and machine learning implementations.

***Primitha Rodrigues*** is an Associate Software Engineer and Data Science Engineer at Optum Technology, specializing in data science and machine learning.

If you found this article helpful, please consider clapping and sharing. Follow us for more deep dives into machine learning fundamentals!