Logistic Regression

Background

Developed in the 1950s-1960s as a way to model binary outcomes in epidemiology and social science
One of the most-used algorithms in industry despite being 70+ years old
Maximum likelihood estimation connects it directly to information theory

The Concept

Logistic regression applies a sigmoid function to a linear combination of features, squashing the output to a probability between 0 and 1. Training maximizes the likelihood of the observed labels, equivalent to minimizing cross-entropy loss.

How it Links

The sigmoid activation in logistic regression is the same one in neural network nodes
Multinomial logistic regression with a softmax output is the final layer of most classifiers
Connects to GLMs, Bayesian inference, and the cross-entropy loss ubiquitous in deep learning

Concrete Example

A bank model predicting loan default: features like income, debt, and payment history get weights; the sigmoid maps their sum to a default probability between 0 and 1. A probability > 0.5 triggers rejection.

Prerequisites

Math you need: linear algebra and probability, Linear regression, Perceptron, Evaluation metrics

Learning outcomes

State the sigmoid model and derive the log-likelihood and its gradient.
Recognize logistic regression as a generalized linear model with the canonical logit link.
Derive Newton-Raphson and IRLS as second-order optimization on the log-likelihood.
Explain why logistic regression is well-calibrated while decision-tree leaf proportions are not.
Implement logistic regression in NumPy (IRLS), sklearn (l1 and l2), and PyTorch (BCEWithLogitsLoss).

TL;DR

Logistic regression models $P(y = 1 \mid x) = \sigma(w^{\top}x + b)$ , where $\sigma$ is the logistic sigmoid. It is the GLM with the canonical logit link and Bernoulli likelihood, and its parameters are fit by maximum likelihood. The negative log-likelihood is convex with gradient $X^{\top}(\hat{p} - y)$ , so it is solved by gradient descent or, more commonly, by Newton-Raphson which becomes iteratively reweighted least squares (IRLS) [1][3]. The decision boundary is still a hyperplane in feature space, but the output is a calibrated probability rather than a hard label as with the perceptron [4].

Intuition

If you want a probability that an example is class 1, logistic regression is the simplest model that gives you one. It has the same affine score $z = w^{\top}x + b$ as the perceptron, but instead of thresholding at zero with the sign function it squashes through the logistic sigmoid

\sigma(z) = \frac{1}{1 + e^{-z}},

which maps $\mathbb{R} \to (0, 1)$ and is monotonic. Two facts make this model well behaved:

The log-odds is linear. The function $\sigma$ inverts to $\sigma^{-1}(p) = \log\frac{p}{1-p}$ , the logit. So logistic regression is exactly the model “the log-odds of class 1 is a linear function of $x$ ” - a small, falsifiable assumption.
The boundary is still a hyperplane. The decision rule $\sigma(z) > 0.5$ is the same as $z > 0$ . The sigmoid does not curve the boundary; it curves the probability surface that sits above it.

Calibration is free here

Logistic regression is trained on the log-likelihood of the observed labels, so its outputs are designed to match the empirical conditional rate $P(y \mid x)$ when the linear-in-log-odds assumption holds. Decision-tree leaf proportions are noisy because each leaf has a small sample, and they are systematically over-confident on the training set. That is why a logistic head on top of a calibrated representation is often the right last layer for production classifiers [1].

Math

Math: log-likelihood, gradient, and Newton step

Model

For binary labels $y \in \{0, 1\}$ , the model is

p(x; w, b) = \sigma(w^{\top}x + b), \qquad \sigma(z) = \frac{1}{1 + e^{-z}}.

A useful identity: $\sigma'(z) = \sigma(z)(1 - \sigma(z))$ .

Likelihood

Assume the labels are conditionally Bernoulli, $y_i \mid x_i \sim \mathrm{Bern}(p(x_i))$ . The likelihood of the whole training set is

\mathcal{L}(w, b) = \prod_{i=1}^{n} p(x_i)^{y_i}\,(1 - p(x_i))^{1 - y_i},

so the log-likelihood is

\ell(w, b) = \sum_{i=1}^{n}\Bigl[\,y_i \log p(x_i) + (1 - y_i)\log(1 - p(x_i))\,\Bigr].

The negative log-likelihood is the binary cross-entropy loss; minimizing it equals maximum likelihood.

Gradient

Let $z_i = w^{\top}x_i + b$ and $p_i = \sigma(z_i)$ . Differentiating term by term, one shows

\frac{\partial \ell}{\partial z_i} = y_i - p_i.

By the chain rule,

\nabla_{w}\,\ell(w, b) = \sum_{i=1}^{n} (y_i - p_i)\, x_i = X^{\top}(y - \hat{p}), \qquad \frac{\partial \ell}{\partial b} = \sum_{i=1}^{n}(y_i - p_i).

The negative log-likelihood gradient is therefore $X^{\top}(\hat{p} - y)$ - the same shape as the OLS gradient, but with the sigmoid squashed predictions in place of the linear ones. There is no closed form, because $\hat{p}$ depends nonlinearly on $w$ .

Hessian

Let $W = \operatorname{diag}\bigl(p_i(1 - p_i)\bigr)$ . Then

\nabla^{2}_{w}\,\bigl[-\ell(w, b)\bigr] = X^{\top} W X \succeq 0,

so the negative log-likelihood is convex (and strictly convex if $X$ has full column rank), guaranteeing a unique minimum and well-behaved optimization.

Newton’s update and IRLS

Newton-Raphson on the negative log-likelihood gives the update

w_{t+1} = w_{t} + (X^{\top} W_t X)^{-1} X^{\top}(y - \hat{p}_t),

which can be rewritten as a weighted least-squares problem:

w_{t+1} = (X^{\top} W_t X)^{-1} X^{\top} W_t z_t, \qquad z_t = X w_t + W_t^{-1}(y - \hat{p}_t).

This is exactly the iteratively reweighted least squares (IRLS) algorithm: each step is a weighted OLS regression of the working response $z_t$ on $X$ with weights $W_t$ . IRLS is the standard fitting procedure for any GLM, not just logistic regression [3].

GLM view

A generalized linear model has three components:

Random component: $y \mid x$ comes from an exponential-family distribution. Here, Bernoulli.
Systematic component: a linear predictor $\eta = w^{\top}x + b$ .
Link function: $g$ relates the mean to the linear predictor via $g(\mu) = \eta$ .

For Bernoulli, the canonical link is the logit, $g(\mu) = \log\frac{\mu}{1 - \mu}$ . Plugging in, $\log\frac{p}{1-p} = w^{\top}x + b$ , which inverts to $p = \sigma(w^{\top}x + b)$ . So logistic regression is the GLM with Bernoulli random component and canonical (logit) link. With the canonical link, the score equation simplifies to $X^{\top}(y - \hat{p}) = 0$ , exactly what we derived [3][2].

Multinomial / softmax extension

For $K$ classes, replace the sigmoid with the softmax

P(y = k \mid x) = \frac{e^{w_{k}^{\top}x + b_{k}}}{\sum_{j=1}^{K} e^{w_{j}^{\top}x + b_{j}}}.

The negative log-likelihood becomes the categorical cross-entropy, and the gradient w.r.t. $w_k$ is $\sum_{i}(\mathbb{1}[y_i = k] - P(y = k \mid x_i))\,x_i$ - the multi-class analogue of $X^{\top}(\hat{p} - y)$ . It is still a convex problem.

Regularized variants

Add an L2 penalty $\frac{\lambda}{2}\lVert w\rVert_{2}^{2}$ for ridge-style logistic regression (MAP under a Gaussian prior on $w$ ); add an L1 penalty $\lambda \lVert w \rVert_{1}$ for sparse logistic regression (MAP under a Laplace prior). Both keep convexity. L2 stabilizes IRLS when $X^{\top}WX$ is ill-conditioned; L1 yields feature selection.

Geometry of the decision boundary

The model predicts class 1 iff $\sigma(w^{\top}x + b) > 0.5$ , which is iff $w^{\top}x + b > 0$ . So the decision boundary is the hyperplane $\{x : w^{\top}x + b = 0\}$ - identical to the perceptron and the linear SVM. The sigmoid is a monotonic post-processing of the score that produces calibrated probabilities; it does not curve the boundary.

What changes between linear classifiers is the loss and therefore which hyperplane is chosen. Perceptron stops at any separating hyperplane; linear SVM picks the max-margin one; logistic regression picks the one that maximises Bernoulli likelihood, which uses all training points, not just the support vectors.

How it compares to nearby linear classifiers

Method	Loss	Output	Calibrated?	Multi-class strategy	Closed form?	Cost per epoch (n examples, d features)
Perceptron	max(0, -y (w·x + b))	class label in {-1, +1}	no	one-vs-rest by hand	no	O(n d)
Logistic regression	binary cross-entropy / negative Bernoulli log-likelihood	P(y = 1 \| x) in (0, 1)	yes (under linear-in-logit assumption)	softmax (multinomial)	no - convex; use IRLS / GD	O(n d) per IRLS iter, O(n d^2 + d^3) for the WLS step
Linear SVM	hinge: max(0, 1 - y (w·x + b))	class label; Platt-scaled probability	no by default; yes after Platt scaling	one-vs-rest or one-vs-one	no - convex QP	O(n d) per SGD epoch, full QP otherwise

Algorithm

# Gradient descent on the negative log-likelihood (binary cross-entropy)
w = zeros(d)
b = 0.0
for step in range(max_steps):
    z = X @ w + b                      # affine score, shape (n,)
    p = 1.0 / (1.0 + exp(-z))          # sigmoid -> probabilities
    grad_w = X.T @ (p - y)             # canonical-link gradient
    grad_b = (p - y).sum()
    w = w - eta * grad_w
    b = b - eta * grad_b
 
# Newton / IRLS: solve a weighted least-squares problem each step
for step in range(max_steps):
    z = X @ w + b
    p = 1.0 / (1.0 + exp(-z))
    W = diag(p * (1 - p))              # working weights
    H = X.T @ W @ X                    # Hessian (positive semidefinite)
    grad = X.T @ (p - y)
    w = w - solve(H, grad)             # Newton step

Implementations

python

import numpy as np

def sigmoid(z):
  # numerically stable sigmoid
  return np.where(z >= 0, 1.0 / (1.0 + np.exp(-z)), np.exp(z) / (1.0 + np.exp(z)))

def logreg_irls(X, y, max_iter=25, tol=1e-8):
  """Fit binary logistic regression by IRLS (Newton on the log-likelihood)."""
  n = X.shape[0]
  X1 = np.hstack([np.ones((n, 1)), X])
  d = X1.shape[1]
  w = np.zeros(d)
  for it in range(max_iter):
      z = X1 @ w
      p = sigmoid(z)
      W = p * (1.0 - p)                                     # diag of W
      # Newton step: w <- w + (X^T W X)^{-1} X^T (y - p)
      H = (X1.T * W) @ X1
      g = X1.T @ (y - p)
      delta = np.linalg.solve(H + 1e-8 * np.eye(d), g)      # tiny ridge for stability
      w_new = w + delta
      if np.linalg.norm(delta) < tol:
          return w_new[0], w_new[1:], it + 1
      w = w_new
  return w[0], w[1:], max_iter

# tiny synthetic dataset
rng = np.random.default_rng(0)
X = rng.normal(size=(200, 2))
y = (X[:, 0] - X[:, 1] + 0.3 * rng.normal(size=200) > 0).astype(float)

b, w, iters = logreg_irls(X, y)
print(f"converged in {iters} Newton steps")
print("intercept :", round(float(b), 3))
print("weights   :", np.round(w, 3).tolist())

IRLS = Newton on the log-likelihood. Note the canonical-link gradient X^T (y - p).

Run the from-scratch NumPy IRLS solver live in the browser:

Idle (packages: numpy)

import numpy as np

def sigmoid(z):
  return np.where(z >= 0, 1.0 / (1.0 + np.exp(-z)), np.exp(z) / (1.0 + np.exp(z)))

def logreg_irls(X, y, max_iter=25, tol=1e-8):
  n = X.shape[0]
  X1 = np.hstack([np.ones((n, 1)), X])
  d = X1.shape[1]
  w = np.zeros(d)
  for it in range(max_iter):
      z = X1 @ w
      p = sigmoid(z)
      W = p * (1.0 - p)
      H = (X1.T * W) @ X1
      g = X1.T @ (y - p)
      delta = np.linalg.solve(H + 1e-8 * np.eye(d), g)
      w = w + delta
      if np.linalg.norm(delta) < tol:
          return w[0], w[1:], it + 1
  return w[0], w[1:], max_iter

rng = np.random.default_rng(0)
X = rng.normal(size=(300, 2))
y = (1.5 * X[:, 0] - X[:, 1] + 0.4 * rng.normal(size=300) > 0).astype(float)

b, w, iters = logreg_irls(X, y)
print("converged in", iters, "Newton steps")
print("intercept :", round(float(b), 3))
print("weights   :", np.round(w, 3).tolist())

# decision-boundary check
preds = (sigmoid(b + X @ w) > 0.5).astype(float)
print("training accuracy :", round(float((preds == y).mean()), 4))

(no output yet)

Real-world uses

Domain	System	Why logistic regression
AdTech / search ranking	Click-through-rate (CTR) baselines	Calibrated probabilities feed downstream auction logic; one-feature-per-coefficient interpretation supports auditability.
Credit / risk	Credit scoring (FICO-style scorecards)	Per-feature coefficients give clear adverse-action explanations required by regulation; closed-form WLS fits at scale.
Healthcare	Clinical prediction models (e.g. Framingham)	Sigmoid output is a probability, not a hand-tuned score; coefficients are interpretable as log-odds-ratios per unit.
NLP	Document classification with bag-of-words / TF-IDF features	Trains in seconds on millions of features; L1-regularized variants double as feature selection.
Calibration head	Last layer of deep classifiers	When the network output is already a logit, BCEWithLogitsLoss is just logistic regression in disguise.

Pitfalls

Perfect separation makes weights blow up

If a feature perfectly separates the two classes, the MLE has $\lVert w \rVert \to \infty$ - the log-likelihood keeps improving without bound. IRLS will fail to converge (the Hessian goes singular). Fix with L2 regularization, Firth’s bias correction, or by replacing logistic with a max-margin loss.

Class imbalance distorts the intercept

The intercept absorbs the marginal class prior. With heavy imbalance the learned threshold of 0.5 is wrong; either use class weights, calibrate the threshold on a validation set, or rescale the intercept by the log prior ratio.

Always standardize before regularizing

L1 and L2 penalties depend on the scale of the features. A feature in cents gets penalized far less than the same feature in dollars. Standardize $X$ before fitting any penalized logistic regression, and never penalize the intercept.

Log-loss is a proper scoring rule

Calibration matters for risk-sensitive applications. Train on log-loss (the binary cross-entropy), then evaluate on log-loss, Brier score, and reliability diagrams - these are all proper scoring rules that reward true probabilities and punish over-confidence.

Interview Q&A

Derive the gradient of the logistic-regression log-likelihood.research

For binary labels $y_i \in \{0, 1\}$ and the model $p_i = \sigma(z_i)$ with $z_i = w^{\top}x_i + b$ , the log-likelihood is

\ell(w, b) = \sum_{i=1}^{n}\bigl[\,y_i \log p_i + (1 - y_i) \log(1 - p_i)\,\bigr].

Use the chain rule. Define $u_i = y_i \log p_i + (1 - y_i)\log(1 - p_i)$ . Then

\frac{\partial u_i}{\partial p_i} = \frac{y_i}{p_i} - \frac{1 - y_i}{1 - p_i} = \frac{y_i - p_i}{p_i (1 - p_i)},

and using $\sigma'(z) = \sigma(z)(1 - \sigma(z))$ ,

\frac{\partial p_i}{\partial z_i} = p_i (1 - p_i).

Multiplying,

\frac{\partial u_i}{\partial z_i} = y_i - p_i,

and finally $\partial z_i / \partial w = x_i$ . Summing,

\nabla_{w}\,\ell(w, b) = \sum_{i=1}^{n} (y_i - p_i)\, x_i = X^{\top}(y - \hat{p}).

The cancellation between $\sigma'(z) = p(1-p)$ in the chain rule and the $1/(p(1-p))$ from $\partial u / \partial p$ is why the gradient is so clean - this is a special property of the canonical-link GLM [3][2].

Why is logistic regression a generalized linear model?onsite

A GLM has three pieces:

Random component: the response distribution is in the exponential family. The Bernoulli is, with natural parameter $\theta = \log\frac{p}{1-p}$ .
Systematic component: a linear predictor $\eta = w^{\top}x + b$ .
Link function $g$ : relates the mean $\mu$ to the linear predictor via $g(\mu) = \eta$ .

For logistic regression, the link is the canonical one for Bernoulli, the logit: $g(\mu) = \log \frac{\mu}{1-\mu}$ . So the model says $\log \frac{p}{1-p} = w^{\top}x + b$ , which inverts to $p = \sigma(w^{\top}x + b)$ . With the canonical link, the score equation is the simple form $X^{\top}(y - \hat{p}) = 0$ , the Hessian is $X^{\top} W X$ with $W$ diagonal, and Newton’s method becomes IRLS - a weighted least-squares problem each step. This is exactly the GLM machinery from McCullagh and Nelder [3].

Why is log-loss a proper scoring rule and what does that mean for training?research

A scoring rule $S(p, y)$ is strictly proper if its expected value over $y \sim q$ is uniquely minimized by $p = q$ . For log-loss,

S(p, y) = -\bigl[y \log p + (1 - y) \log (1 - p)\bigr],

we compute the expected score under the true distribution $q$ :

\mathbb{E}_{y \sim q}\,S(p, y) = -\bigl[q \log p + (1 - q)\log(1 - p)\bigr] = H(q) + \mathrm{KL}(q\,\Vert\,p),

which is minimized at $p = q$ because $\mathrm{KL}(q\,\Vert\,p) \ge 0$ with equality iff $p = q$ . So if the model class is rich enough to represent the true conditional $P(y \mid x)$ , fitting on log-loss recovers it - log-loss penalizes mis-calibrated probabilities, not just mis-classifications. This is the foundation of why logistic regression outputs are calibrated probabilities. It is also why people use log-loss to compare classifiers when calibration matters: a 0/1-accuracy comparison can rank a worse-calibrated model higher [1].

Derive Newton's update on the negative log-likelihood and show it equals IRLS.research

The negative log-likelihood is $\mathrm{NLL}(w) = -\ell(w)$ . We computed

\nabla_{w}\,\mathrm{NLL}(w) = -X^{\top}(y - \hat{p}) = X^{\top}(\hat{p} - y),

\nabla^{2}_{w}\,\mathrm{NLL}(w) = X^{\top} W X, \qquad W = \operatorname{diag}\bigl(p_i(1 - p_i)\bigr).

Newton’s step is

w_{t+1} = w_{t} - \bigl(X^{\top} W_t X\bigr)^{-1} X^{\top}(\hat{p}_t - y) = w_{t} + \bigl(X^{\top} W_t X\bigr)^{-1} X^{\top}(y - \hat{p}_t).

Define the working response $z_t = X w_t + W_t^{-1}(y - \hat{p}_t)$ (the linearization of the link). Plugging in,

w_{t+1} = \bigl(X^{\top} W_t X\bigr)^{-1} X^{\top} W_t z_t,

which is the closed-form solution to the weighted least-squares problem $\min_{w} (z_t - X w)^{\top} W_t (z_t - X w)$ . Each Newton step is therefore a weighted OLS regression of the working response on the features, with weights $W_t$ updated each iteration. That is iteratively reweighted least squares, which is the fitting algorithm for GLMs and is what glm() in R and statsmodels.GLM actually run under the hood [3].

Check yourself

Score: 0 / 5(0 answered)

1. What is the gradient of the negative log-likelihood for binary logistic regression?
2. Why does the negative log-likelihood for logistic regression have NO closed form?
3. In the GLM framework, logistic regression uses which canonical link?
4. What is the geometry of the decision boundary for logistic regression?
5. Why are logistic-regression probabilities considered calibrated, while decision-tree leaf proportions are not?

Build it yourself

A short evening project:

Generate a 2-class synthetic dataset with overlap (e.g. two Gaussians with overlapping means).
Fit logistic regression by IRLS, by gradient descent, and by sklearn; confirm all three give the same coefficients to within 1e-6.
Plot the decision boundary and the probability heat map.
Add 50% L2 regularization and watch the boundary tilt back towards the origin.
Stretch: replicate Platt scaling - take a hard-margin linear SVM, then fit a logistic regression on its scores to recover calibrated probabilities, and compare to native logistic regression on log-loss.

Interactive Demo

Logistic Regression Boundary Demo

w₁ 1.00w₂ -1.00b 0.00

P(y=1) = σ(w₁x + w₂y + b)

Cross-entropy loss: 0.0528

Class A (red)

Class B (blue)

Decision boundary

What to try: Adjust the weight sliders to see how the decision boundary and probability heatmap shift.

References

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed.. Springer. https://hastie.su.domains/ElemStatLearn/
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2
McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models, 2nd ed.. Chapman & Hall / CRC. ISBN 978-0-412-31760-6
Rosenblatt, F. (1958). The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain. Psychological Review 65(6), 386-408. doi:10.1037/h0042519

Linear Regression k-Nearest Neighbors