Logistic Regression from Scratch in Python: A Deep Dive

Logistic regression is a foundational machine learning algorithm for binary classification problems. While it‘s easy to use powerful libraries like scikit-learn to train logistic regression models, building the algorithm from scratch is one of the best ways to understand how it really works under the hood.

In this in-depth guide, we‘ll derive logistic regression from first principles and implement it in Python, covering the underlying math, optimization process, and extensions in detail. Mastering these fundamentals will give you a strong intuition for training and debugging machine learning models.

Deriving Logistic Regression

Let‘s start by deriving the logistic regression model mathematically. Assume we have input features X = [x1, x2, …, xm] and binary labels y ∈ {0, 1}. We want to model the probability that y=1 given X.

A natural choice is the Bernoulli distribution, which models a single binary outcome with probability p:

P(y|p) = p^y * (1-p)^(1-y)

We can parameterize p as a function of the inputs X using the sigmoid function σ(z) = 1/(1+e^(-z)):

p = σ(θ^T X) = 1/(1+e^(-θ^T X))

Here, θ are the model‘s weight parameters. The sigmoid "squashes" θ^T X to a probability between 0 and 1.

Plugging this into the Bernoulli distribution, we get the likelihood of the data under the model:

L(θ) = ∏ σ(θ^T Xi)^yi * (1-σ(θ^T Xi))^(1-yi)

To simplify optimization, we take the log likelihood:

ℓℓ(θ) = ∑ yi·log(σ(θ^T Xi)) + (1-yi)log(1-σ(θ^T Xi))

We want to find the parameters θ that maximize the log likelihood, which is equivalent to minimizing the negative log likelihood (NLL). This is our loss function J(θ):

J(θ) = -1/m ∑ yi·log(σ(θ^T Xi)) + (1-yi)log(1-σ(θ^T Xi))

Minimizing J(θ) will give us the optimal parameters θ that fit the data. This loss is also called the cross-entropy or log loss.

Optimization via Gradient Descent

To minimize the loss J(θ), we‘ll use batch gradient descent. The gradient ∇J(θ) tells us the direction of steepest ascent. By taking steps in the negative gradient direction, we can iteratively descend the loss curve toward a minimum.

The jth component of the gradient is the partial derivative of J(θ) with respect to θj:

∂J(θ)/∂θj = -1/m ∑ (yi – σ(θ^T Xi)) · Xij

We can derive this using the chain rule. The full gradient ∇J(θ) is then:

∇J(θ) = -1/m X^T · (y – σ(Xθ))

Where y and σ(Xθ) are m-dimensional column vectors.

The gradient descent update rule is:

θ := θ – α · ∇J(θ)

Where α is the learning rate that controls the step size.

Intuitively, if the model predicts a probability σ(θ^T Xi) that‘s very different from the true label yi, the gradient update will be large, causing θ to change significantly. Conversely, if the prediction is already good, the update will be small.

Here‘s a visualization of gradient descent on the logistic loss curve for a single parameter θ:

The blue curve is the log loss as a function of θ. The red lines show the negative gradient -∇J(θ) at different points. Gradient descent iteratively takes steps in this direction, moving θ toward the global minimum.

Python Implementation

Let‘s put the pieces together and implement logistic regression in Python.

We‘ll start by importing NumPy and generating a synthetic dataset using scikit-learn‘s make_classification function:

import numpy as np
from sklearn.datasets import make_classification

X, y = make_classification(n_samples=1000, n_features=5, random_state=42)

Next, we‘ll define the sigmoid function, loss function, and gradient function:

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

def loss(X, y, w):
    z = X @ w
    p = sigmoid(z)
    return -np.mean(y*np.log(p) + (1-y)*np.log(1-p))

def gradient(X, y, w):
    z = X @ w
    p = sigmoid(z)
    return (1/len(y)) * X.T @ (p - y)

We can now implement gradient descent:

def logistic_regression(X, y, learning_rate=0.1, num_iterations=10000):
    w = np.zeros(X.shape[1])
    loss_history = []

    for i in range(num_iterations):
        w -= learning_rate * gradient(X, y, w)
        loss_history.append(loss(X, y, w))

    return w, loss_history

Finally, we‘ll train the model and plot the training loss over time:

w, loss_history = logistic_regression(X, y)

plt.plot(loss_history)
plt.xlabel(‘Iteration‘)
plt.ylabel(‘Log Loss‘)
plt.title(‘Training Loss‘)
plt.show()

The loss decreases rapidly at first and then plateaus as the model converges.

We can make predictions on new data using our learned weights:

def predict(X, w):
    return np.round(sigmoid(X @ w))

y_pred = predict(X, w)
print(f‘Training accuracy: {np.mean(y_pred == y):.3f}‘)

Training accuracy: 0.901

Our from-scratch model achieves over 90% accuracy on the synthetic data, which is great considering we haven‘t done any feature engineering or hyperparameter tuning.

Regularization

One way to improve logistic regression is to add regularization, which penalizes large weight values to mitigate overfitting. The two common types are L1 (Lasso) and L2 (Ridge) regularization.

L1 regularization adds the absolute values of the weights to the loss:

J(θ) = -1/m ∑ yi·log(σ(θ^T Xi)) + (1-yi)log(1-σ(θ^T Xi)) + λ∑|θj|

While L2 regularization adds the squared values:

J(θ) = -1/m ∑ yi·log(σ(θ^T Xi)) + (1-yi)log(1-σ(θ^T Xi)) + λ∑θj^2

Where λ controls the regularization strength. Larger λ values apply stronger penalties, encouraging the model to learn smaller weights.

We can implement L2 regularization in our loss and gradient functions like this:

def loss(X, y, w, lam=1.0):
    z = X @ w
    p = sigmoid(z)
    return -np.mean(y*np.log(p) + (1-y)*np.log(1-p)) + (lam/(2*len(y))) * np.sum(w**2)

def gradient(X, y, w, lam=1.0):
    z = X @ w
    p = sigmoid(z)
    return (1/len(y)) * X.T @ (p - y) + (lam/len(y)) * w

We can visualize the impact of different λ values on the learned weights using a regularization path:

lambdas = np.logspace(-3, 3, 100)
weights = []

for lam in lambdas:
    w, _ = logistic_regression(X, y, learning_rate=0.1, num_iterations=10000, lam=lam)
    weights.append(w)

plt.plot(lambdas, weights)
plt.xscale(‘log‘)
plt.xlabel(‘Lambda‘)
plt.ylabel(‘Weights‘)
plt.title(‘Regularization Path‘)
plt.show()

As λ increases, the weights shrink toward zero. Regularization is a powerful tool for controlling model complexity and preventing overfitting.

Comparative Benchmarks

Let‘s compare our custom implementation to scikit-learn‘s on a few standard datasets:

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

datasets = [make_classification(n_samples=1000, n_features=20, random_state=42),
            make_classification(n_samples=10000, n_features=100, random_state=42),
            make_classification(n_samples=100000, n_features=50, n_informative=10, random_state=42)]

for X, y in datasets:
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

    # Custom model
    w, _ = logistic_regression(X_train, y_train)
    y_pred = predict(X_test, w)
    print(f‘Custom model accuracy: {accuracy_score(y_test, y_pred):.3f}‘)

    # sklearn model  
    clf = LogisticRegression(penalty=‘none‘)
    clf.fit(X_train, y_train)
    print(f‘sklearn model accuracy: {clf.score(X_test, y_test):.3f}\n‘)

Custom model accuracy: 0.895
sklearn model accuracy: 0.895

Custom model accuracy: 0.888 
sklearn model accuracy: 0.888

Custom model accuracy: 0.733
sklearn model accuracy: 0.734

Our custom implementation achieves near-identical performance to scikit-learn‘s on all three datasets, which is impressive! However, scikit-learn‘s version is much more computationally efficient, leveraging compiled C code under the hood.

In practice, it‘s recommended to use optimized library implementations for production use cases. However, implementing algorithms from scratch is invaluable for deepening understanding.

Multiclass Logistic Regression

So far, we‘ve focused on binary classification. To extend logistic regression to multiclass problems, we use the "one-vs-rest" (OvR) approach, also known as "one-hot encoding".

The idea is to train K independent binary classifiers, one for each of the K classes. The kth classifier is trained to predict the probability that an instance belongs to class k vs. not-k.

Concretely, for a 3-class problem with labels [cat, dog, bird], we‘d train 3 binary classifiers:

1. P(y = cat | X)
2. P(y = dog | X)
3. P(y = bird | X)

To make a prediction, we run an instance through all K classifiers and choose the class with the highest probability.

In practice, we can efficiently implement OvR multiclass logistic regression using the softmax function. Softmax extends the sigmoid to the multiclass case:

P(y=k|X) = e^(θ_k^T X) / ∑e^(θ_j^T X)

Where θ_k are the weights for the kth class. The softmax normalizes the exponentiated scores into a valid probability distribution over the K classes.

Most ML libraries like scikit-learn handle this under the hood, but it‘s good to understand the underlying mechanics.

Extensions & Alternatives

Logistic regression is a special case of a generalized linear model (GLM) with the logit link function. Other GLMs like Poisson regression (for count data) and Gamma regression (for strictly positive data) use different link functions and loss distributions.

An alternative to the sigmoid is the probit function Φ(z), which uses the CDF of the standard normal distribution instead of the logistic function. In practice, the probit and logit are very similar, with the probit having slightly heavier tails.

More advanced extensions include:

• Elastic net regularization: Combines L1 and L2 penalties for sparse and stable solutions
• Bayesian logistic regression: Places priors on weights to incorporate domain knowledge
• Kernel logistic regression: Applies the kernel trick for non-linear decision boundaries
• Online learning: Stochastic gradient descent for incremental learning on streaming data

Diving into these extensions is a great way to build on your understanding of core logistic regression.

Conclusion

We covered a lot of ground in this deep dive! Starting from the fundamentals of binary classification and maximum likelihood estimation, we derived the logistic regression model and optimization objective.

We then implemented the core components of logistic regression from scratch in Python, including batch gradient descent, L2 regularization, and comparative benchmarks to scikit-learn.

Key takeaways:

• Logistic regression is a discriminative classifier that models P(y|X) using the Bernoulli distribution and sigmoid link function
• Log loss is the principled probabilistic loss function for binary data and enables gradient-based optimization
• Regularization mitigates overfitting by penalizing large weights
• Multiclass logistic regression uses the softmax function and one-vs-rest strategy
• Logistic regression is a special case of the broader family of generalized linear models

I hope this guide clarified the inner workings of logistic regression and equipped you with a solid foundation for further exploration. Implementing algorithms from scratch is a powerful learning tool that I highly recommend for all fundamental ML models.

The complete code is available on GitHub. Feel free to use it as a starting point for your own projects and experiments. If you have any questions or insights to share, please leave a comment below. Happy learning!

References

• Bishop, C. (2006). Pattern Recognition and Machine Learning. Springer.
• Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning. MIT Press.
• Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.). Springer.
• Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press.
• Ng, A. (2012). CS229 Lecture Notes. Stanford University.