how to implement logistic regression from scratch

What is logistic regression?

Welcome back! Today we are tackling a concept that sounds intimidating but is actually quite elegant. Logistic regression is a workhorse algorithm for binary classification.

Despite the name, it doesn't predict continuous numbers like Linear Regression does. Instead, it predicts the probability that an input belongs to the positive class (often labeled as 1). The output is a number strictly between 0 and 1. We can then take this probability and apply a threshold (usually 0.5) to make a final decision.

Misconception: Is it only for linear data?

It is true that logistic regression assumes a linear relationship between the features and the log-odds (the logarithm of the odds). However, this doesn't mean your data must be perfectly linear.

By engineering new features—like adding polynomial terms (e.g., $x^2$) or interactions ($x_1 \times x_2$)—you can capture complex, non-linear patterns while still using the logistic regression machinery. The model remains linear in its parameters, but flexible enough for complex data.

How the S-shaped curve emerges

The magic happens through the Sigmoid function (also called the logistic function). This function takes any real number (from $-\infty$ to $+\infty$) and "squashes" it into a range between 0 and 1.

Here is how you would implement this squashing function in Python:

import numpy as np

def sigmoid(z):
    # The magic formula that maps any number to 0-1
    return 1 / (1 + np.exp(-z))

# Try it out!
print(sigmoid(0))    # Output: 0.5
print(sigmoid(10))   # Output: ~0.999
print(sigmoid(-10))  # Output: ~0.000

This S-curve is the core of logistic regression. It maps the raw score (linear combination of weights and inputs) to a probability, enabling us to make binary classification decisions.

Binary classification with logistic regression

Now that we understand the Sigmoid curve, let's look at how logistic regression actually uses it to solve problems. Imagine you have a scatter plot with two groups of points: Blue Circles (Class 0) and Red Crosses (Class 1).

Misconception: Does it output "Yes" or "No"?

This is a very common point of confusion. Logistic regression does not output a class label like "Spam" or "Not Spam" directly.

Its raw output is always a probability ($p$) between 0 and 1. It is you (the developer) who decides how to interpret that number. You choose a threshold (usually 0.5) to convert the probability into a final class label.

def predict_labels(probabilities, threshold=0.5):
    # Convert probability to binary label
    return (probabilities >= threshold).astype(int)

# Example usage
prob = 0.75
label = predict_labels(prob) 
# Output: 1 (since 0.75 >= 0.5)

The default threshold of 0.5 works well when both classes are equally important. However, in the real world, you might need to adjust this "knob". For example, in cancer detection, missing a positive case (False Negative) is very dangerous. You might lower the threshold to 0.3 so the model predicts "Cancer" more aggressively, ensuring you catch every possible case, even if it means some false alarms.

Why implement machine learning from scratch?

We have discussed the math and the intuition. Now, a crucial question arises: Why bother writing the code yourself? Why not just import a library like scikit-learn or TensorFlow and be done with it?

The short answer is: Control and Clarity. When you implement an algorithm from scratch, you stop treating the model as a magical "black box." You begin to see the gears turning inside.

Misconception: From-scratch implementations are slower

It is absolutely true that libraries like scikit-learn are heavily optimized, often using C++ or CUDA under the hood. They will almost always be faster than a pure Python implementation you write in a notebook.

However, the goal of a from-scratch implementation is not raw speed—it is comprehension. For learning purposes and small datasets, the performance difference is negligible. What you gain is the ability to modify any part of the algorithm.

So, implementing from scratch isn't about replacing professional tools. It's about mastering them. Once you know how the engine works, you can drive the car faster and fix it better when it stalls.

Setting up Python for logistic regression

Welcome to the workshop! Before we start building our algorithm, we need to gather our tools. You are about to build the entire logistic regression logic yourself—no magic, no hidden layers.

The last thing you want is to get tangled in complex installation procedures or heavyweight frameworks before you even start coding. A minimal setup with just a few core scientific libraries keeps your focus on the algorithm's logic, not on configuration headaches.

Installing the essentials

Many beginners assume that machine learning requires TensorFlow or PyTorch right away. That is not true for logistic regression. Those frameworks shine for large-scale deep learning, but here we only need efficient array operations and math functions.

Open your terminal or command prompt and run the following command. That's it!

# Install the core scientific libraries
pip install numpy scipy

Once installed, you will import them in your script. We use NumPy for all array manipulations and element-wise math because it's fast and expressive.

import numpy as np
from scipy import optimize  # We'll use this for gradient descent later

You can verify the installation by checking the version. This ensures everything is connected correctly before we start coding the logic.

print(np.__version__)  # e.g., 1.26.0
print(optimize.__doc__[:100])  # prints a short description

Starting with these lightweight tools ensures you see exactly what each operation does, building a solid foundation. Later, if you choose, you can transition to bigger libraries, but for now, let's keep it simple and transparent.

Step-by-step implementation of the logistic regression algorithm

Welcome to the workshop! Now that we have our tools, let's start building. Implementing an algorithm from scratch is like constructing a house. You don't just throw bricks on the ground; you follow a strict sequence:

Misconception: You can skip data normalization

It is tempting to think that because logistic regression is simple, we can feed it raw data directly. While the math technically works, it is a terrible idea for optimization.

If one feature ranges from 0 to 1 (like probability) and another ranges from 0 to 1000 (like income), the "Cost Function" landscape becomes a long, narrow valley. Gradient descent will get stuck zig-zagging down the sides of this valley instead of going straight to the bottom.

Data preprocessing steps

Let's write the code. We assume you have a CSV file. We will use NumPy for the heavy lifting and sklearn just for the split (though we could do that manually too).

import numpy as np

# Load data (skip header row)
data = np.loadtxt('data.csv', delimiter=',', skiprows=1)

# Separate features (X) and labels (y)
X = data[:, :-1]  # All columns except last
y = data[:, -1]   # Last column

# Create a mask for rows without NaNs
mask = ~np.isnan(X).any(axis=1)

# Filter X and y
X = X[mask]
y = y[mask]

# Calculate mean and std
mean = X.mean(axis=0)
std = X.std(axis=0)

# Apply standardization
X_scaled = (X - mean) / std

# Append a column of ones to the start
X_b = np.c_[np.ones((X_scaled.shape[0], 1)), X_scaled]
# Shape is now (m, n+1)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X_b, y, test_size=0.2, random_state=42
)

Once you have X_train, you are ready to define your Sigmoid function and start the training loop. Remember: in a real project, you would fit the scaler on the training set, then use those same statistics to transform the test set.

Gradient descent intuition for logistic regression (Advanced)

Now we arrive at the heart of the matter: how the model actually learns. We have the math, and we have the data. But how do we find the specific numbers (weights) that make the predictions accurate?

We use an algorithm called Gradient Descent. To understand it, we need to visualize the "Cost Function" $J(w)$ not as a formula, but as a physical landscape.

Misconception: Gradient descent always converges quickly and reliably

It is a common belief that because Logistic Regression has a convex cost function (a perfect bowl shape), Gradient Descent will always find the solution instantly. While true that it will find the global minimum eventually, the speed and stability depend entirely on your Learning Rate.

If the learning rate is too large, you take giant leaps. You might jump from one side of the valley to the other, never settling at the bottom. If it is too small, you might take a million tiny steps and never finish your training before the deadline.

Code: Gradient computation and weight update

Let's look at the actual math behind the movement. The gradient of the log loss with respect to weights $w$ is:

Here is a vectorized implementation in Python using NumPy. This code performs the "step" we visualized in the chart above.

import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def compute_gradient(X, y, w):
    """
    X: (m, n+1) augmented feature matrix
    y: (m,) binary labels {0,1}
    w: (n+1,) weight vector
    Returns: (n+1,) gradient vector
    """
    m = X.shape[0]
    z = X.dot(w)          # Linear combination: (m,)
    p = sigmoid(z)        # Predictions: (m,)
    
    # The gradient is the average error weighted by features
    gradient = (1/m) * X.T.dot(p - y)  
    return gradient

def train_logistic_regression(X, y, lr=0.01, epochs=1000, tol=1e-6):
    """
    X: (m, n+1) augmented training data
    y: (m,) training labels
    lr: learning rate
    epochs: max iterations
    tol: convergence tolerance on weight change
    Returns: learned weights w (n+1,)
    """
    n = X.shape[1]
    w = np.zeros(n)  # Initialize weights to zero
    
    for i in range(epochs):
        grad = compute_gradient(X, y, w)
        w_new = w - lr * grad  # Gradient descent step
        
        # Check convergence: if weights barely changed, stop early
        if np.linalg.norm(w_new - w) < tol:
            print(f"Converged at epoch {i}")
            break
        w = w_new
        
    return w

Key points in the code:

• X includes a column of ones for the bias term; w[0] is the intercept.
• The gradient formula is derived by differentiating the log loss; the term p - y is the residual error for each example.
• The update w = w - lr * grad moves opposite the gradient (downhill).
• We check convergence by measuring the Euclidean norm of the weight change; if it's tiny, we've reached the minimum.

With this loop, you now have a complete training procedure. You are essentially programming the model to "feel" the error and correct itself iteratively. In the next section, we'll see how to use these learned weights to make predictions on new data.

Evaluating logistic regression model performance

Welcome back! You have trained your model, and it's time for the report card. But in the world of machine learning, getting a "good grade" is more nuanced than just a single number.

After training your logistic regression model, you need to know how well it performs on unseen data. For binary classification, the most common metrics are accuracy, precision, and recall.

Misconception: High accuracy means a good model

It is easy to celebrate a high accuracy score, but this can be dangerously misleading. Accuracy is not a reliable metric when your classes are imbalanced.

Suppose 95% of emails are not spam (negative class). A naive model that always predicts "not spam" will achieve 95% accuracy, yet it is completely useless—it never detects spam.

Confusion matrix basics

To see the full picture, we use the Confusion Matrix. It is a simple 2x2 table that compares true labels against predicted labels. It counts four specific outcomes:

From these four numbers, you compute all other metrics. Let's play with the numbers to see how they affect the scores.

Code implementation

Now, let's implement these calculations in Python. We assume you have true binary labels y_true and predicted binary labels y_pred (0 or 1) after thresholding probabilities.

import numpy as np

def confusion_matrix(y_true, y_pred):
    """
    Returns counts: TN, FP, FN, TP in that order.
    """
    tn = np.sum((y_true == 0) & (y_pred == 0))
    fp = np.sum((y_true == 0) & (y_pred == 1))
    fn = np.sum((y_true == 1) & (y_pred == 0))
    tp = np.sum((y_true == 1) & (y_pred == 1))
    return tn, fp, fn, tp

def evaluate_model(y_true, y_pred):
    tn, fp, fn, tp = confusion_matrix(y_true, y_pred)
    
    # Calculate metrics with division-by-zero protection
    total = tp + tn + fp + fn
    accuracy = (tp + tn) / total if total > 0 else 0
    
    precision = tp / (tp + fp) if (tp + fp) > 0 else 0
    recall = tp / (tp + fn) if (tp + fn) > 0 else 0
    
    return {
        'accuracy': accuracy,
        'precision': precision,
        'recall': recall,
        'confusion': {'TN': tn, 'FP': fp, 'FN': fn, 'TP': tp}
    }

Why this code matters:

• We compute the four counts using boolean masks on NumPy arrays—this is fast and vectorized.
• Each metric includes a guard against division by zero (e.g., if there are no predicted positives, precision is set to 0). This prevents runtime errors and gives a sensible interpretation.
• The function returns a dictionary with all key metrics and the raw confusion counts, so you can inspect exactly where errors occur.

Next steps: Interpreting the results

Once you have these numbers, what do you do?

Remember, no single metric tells the whole story. Always look at the full picture of the confusion matrix to diagnose your model's behavior.

Common pitfalls and debugging tips

Welcome back to the workshop! You have the math, the code, and the intuition. But in the real world of software engineering, code rarely works perfectly on the first try.

When implementing algorithms from scratch, you often run into "silent killers"—bugs that don't crash your program but produce garbage results. Today, we are going to learn how to spot them.

The Overflow Cliff

Let's visualize why extreme values break the math. The chart below shows the relationship between the input $z$, the sigmoid probability, and the resulting Loss (Error).

Misconception: The model will always converge

It is a common belief that because logistic regression is "convex" (a simple bowl shape), Gradient Descent will always find the solution.

While mathematically true that a solution exists, computationally, things can go wrong. If your learning rate is too high, you might bounce around the bottom of the bowl forever. If your data is perfectly separable, the weights might grow to infinity in an attempt to achieve 100% accuracy.

Code: The "Numerically Stable" Sigmoid

Instead of just clipping the output, we can write the sigmoid function itself to be stable. This version handles large positive and negative numbers safely.

def stable_sigmoid(z):
    """
    Numerically stable sigmoid function.
    Handles large positive and negative z values without overflow.
    """
    # Split into positive and negative cases
    pos_mask = z >= 0
    neg_mask = z < 0
    
    # Initialize result array
    res = np.zeros_like(z, dtype=float)
    
    # For positive z: use standard formula
    # exp(-z) is safe for positive z
    res[pos_mask] = 1.0 / (1.0 + np.exp(-z[pos_mask]))
    
    # For negative z: use equivalent formula to avoid overflow
    # exp(z) is safe for negative z (it becomes small)
    # 1 / (1 + e^-z) = e^z / (1 + e^z)
    exp_z = np.exp(z[neg_mask])
    res[neg_mask] = exp_z / (1.0 + exp_z)
    
    return res

By implementing these checks and using stable math, you transform yourself from a "coder" into a "debugger". You are now ready to handle real-world data, where things are messy and numbers are extreme.

Extending to multi-class (softmax) – (Advanced)

Welcome back! So far, we have been binary thinkers: Yes or No, 0 or 1. But the real world is rarely that simple. What if you need to classify an image as a Cat, a Dog, or a Bird?

To solve this, we generalize our approach. We replace the single Sigmoid function with the Softmax function. Think of Softmax not as a single switch, but as a distributing mechanism that takes raw scores for multiple classes and turns them into a probability distribution.

Code: The Numerically Stable Softmax

Implementing Softmax in Python requires one crucial trick: numerical stability. If you compute $e^{1000}$, the number is too large for a computer to handle (overflow). We fix this by subtracting the maximum value from all scores. This doesn't change the math, but it keeps the numbers safe.

import numpy as np

def softmax(z):
    """
    z: Input array of shape (m, C) where m is examples and C is classes.
    """
    # Subtract max for numerical stability (prevents overflow)
    # This is safe because e^(z - max(z)) = e^z / e^max(z)
    z_max = np.max(z, axis=1, keepdims=True)
    exp_z = np.exp(z - z_max)
    
    # Normalize
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

# Example usage:
# Scores for 3 classes
scores = np.array([[1.0, 2.0, 3.0]]) 
probs = softmax(scores)
# Output: [0.09, 0.24, 0.67] -> Sum is 1.0

Misconception: Logistic Regression is only binary

The name "Logistic Regression" usually implies the binary case (Sigmoid). When we use Softmax, we are technically performing Multinomial Logistic Regression.

Some beginners try to use One-vs-Rest (training one binary classifier per class). While this works, it is different from Softmax. One-vs-Rest treats classes independently. Softmax treats them as a competition.

You have now mastered the transition from binary to multi-class thinking. You understand that Softmax is just a "normalized vote" and that the math behind it is an extension of the same principles we learned with Logistic Regression.

Frequently Asked Questions

You have built the model, trained it, and evaluated it. But you likely have questions. Here are the most common hurdles students face when implementing Logistic Regression, along with the "Professor's" perspective on how to solve them.