Logistic Regression
The sigmoid, log-odds, cross-entropy loss, and the Hessian that proves it's convex.
On this page
Beginner: OLS (section 2.3.1) predicts a number that can be anything from −∞ to +∞. That's fine for house prices; it's useless for "is this email spam" because a probability has to live between 0 and 1. Logistic regression's fix is simple to state: compute the same linear score z = θᵀx as before, then squash it through an S-shaped curve — the sigmoid — that maps any real number into (0, 1). Large positive scores get squashed close to 1, large negative scores close to 0, and a score of exactly 0 lands right on 0.5. Fitting the model then means choosing θ so that this squashed score is close to 1 for the "yes" examples and close to 0 for the "no" examples.
Intermediate: "close to" needs a loss function, and the obvious first guess — squared error between the sigmoid output and the 0/1 label, exactly like OLS — turns out to be the wrong choice. The sigmoid saturates: far from the decision boundary its slope is nearly flat, so squared error's gradient (which carries a stray factor of that slope) goes nearly to zero exactly on the confidently-wrong predictions that most need correcting. Binary cross-entropy, derived below directly from treating the label as a Bernoulli random variable and maximizing its likelihood (section 2.1.6), doesn't have this problem — its gradient collapses that saturating slope term away entirely, leaving an error signal that stays proportional to how wrong the prediction is, no matter how confident.
Advanced: logistic regression already made a brief appearance in section 2.3.9 as the canonical logit-link Bernoulli GLM, fit there by IRLS. This lesson is the same model derived the other way around — bottom-up, from the sigmoid and the Bernoulli likelihood directly to a hand-derived gradient and Hessian — rather than top-down from the GLM machinery. The payoff of doing it this way is the Hessian itself: H = XᵀWX is positive semidefinite everywhere, which by section 2.2.1's characterization of convexity means this loss has no local minima to get stuck in at all — a guarantee gradient descent almost never gets for free, and one that softmax regression (section 2.4.2) inherits directly when this exact derivation is generalized to K classes.
The linear score, unchanged from OLS:
squashed into a probability by the sigmoid:
fit by minimizing the binary cross-entropy (negative log-likelihood):
Every step from "linear score" to "this exact loss" is derived from scratch below — nothing here is asserted without a reason.
Start from the definition, σ(z) = 1/(1+e^{-z}). Three properties make it the right squashing function for this job.
Range. As z → −∞, e^{-z} → ∞, so σ(z) → 0. As z → +∞, e^{-z} → 0, so σ(z) → 1. At z = 0, σ(0) = 1/2. So σ maps all of ℝ onto the open interval (0, 1) — exactly what a probability needs.
Symmetry. σ(-z) = 1/(1+e^{z}). Multiply numerator and denominator by e^{-z}: σ(-z) = e^{-z}/(e^{-z}+1) = 1 - 1/(1+e^{-z}) = 1 - σ(z). So flipping the sign of the score exactly flips the predicted probability across 0.5 — σ(-z) = 1 - σ(z), which is why "predict class 0" and "predict class 1" are perfectly symmetric operations under this model.
Derivative. Write σ(z) = (1+e^{-z})^{-1} and apply the chain rule:
Rewrite the numerator as e^{-z} = (1+e^{-z}) - 1:
This identity, σ'(z) = σ(z)(1-σ(z)), is the single fact both derivations below lean on — it lets every later derivative be written purely in terms of σ itself, with no leftover exponentials.
The log-odds. Since σ is strictly increasing, it has an inverse — solve p = σ(z) for z. From p = 1/(1+e^{-z}): 1/p = 1+e^{-z}, so e^{-z} = (1-p)/p, so z = log(p/(1-p)). This quantity, the log of the odds p/(1-p), is called the logit:
Where this is used: this is exactly the "canonical link function" from section 2.3.9's GLM framing — logistic regression doesn't model P(y=1|x) as linear in x directly (that could produce probabilities outside [0,1]); it models the log-odds as linear in x, and the sigmoid is nothing more than that statement solved back for p.
Treat each label yᵢ ∈ {0, 1} as one draw from a Bernoulli random variable whose success probability the model supplies: pᵢ = σ(θᵀxᵢ). A single compact expression covers both cases at once:
Check it: when yᵢ = 1 this reduces to pᵢ; when yᵢ = 0 it reduces to 1-pᵢ — exactly the two Bernoulli outcomes, written as one formula so it can be differentiated without a case split.
Assuming the n examples are drawn independently, the full likelihood is the product of these terms, and — following the MLE principle from section 2.1.6 exactly — fitting θ means maximizing this likelihood, which is equivalent to maximizing its log (a strictly increasing function doesn't change where the maximum sits):
Optimizers are conventionally written as minimization problems, so define the loss as the negative of this: L(θ) = −log 𝓛(θ). That single sign flip turns "maximize the log-likelihood" into "minimize the negative log-likelihood," and the result is precisely the binary cross-entropy from the Formula section above:
Where this is used: this is why squared error was the wrong first guess in the plain-English section above — it isn't the negative log-likelihood of any Bernoulli model, so minimizing it doesn't correspond to maximum likelihood estimation at all. Cross-entropy is the loss MLE actually produces for this exact likelihood, and (as the next derivation shows directly) its gradient is exactly what fixes squared error's saturating-gradient problem.
Work with one example's contribution to the loss first, zᵢ = θᵀxᵢ, lᵢ(θ) = −yᵢ log σ(zᵢ) − (1−yᵢ) log(1−σ(zᵢ)), and use the chain rule through zᵢ:
since zᵢ = θᵀxᵢ is linear in θ with slope xᵢ. Differentiate lᵢ with respect to zᵢ using the identity σ'(z) = σ(z)(1-σ(z)) derived above:
Expand and collect terms: -y_i + y_i\sigma(z_i) + \sigma(z_i) - y_i\sigma(z_i) — the two y_i\sigma(z_i) terms cancel exactly, leaving:
That cancellation is the whole reason cross-entropy fixes squared error's saturation problem: the sigmoid's own slope, σ', disappears completely from the gradient, leaving a bare prediction-error term (pᵢ − yᵢ) that stays large exactly when the model is confidently wrong. Substituting back:
Summing over all n examples and writing the result in matrix form:
— the clean closed form the content scope for this lesson promises: stack the per-example errors σ(Xθ) − y into one vector and left-multiply by Xᵀ, exactly mirroring OLS's normal-equations gradient Xᵀ(Xθ − y) from section 2.3.1, but with the linear prediction Xθ replaced by the squashed prediction σ(Xθ).
Now the Hessian. Differentiate ∇L once more. Each term (σ(θᵀxᵢ) − yᵢ)xᵢ depends on θ only through σ(θᵀxᵢ), and by the chain rule again, \partial \sigma(\theta^\top x_i)/\partial \theta = \sigma'(\theta^\top x_i)\,x_i = \sigma_i(1-\sigma_i)\,x_i. So the derivative of the i-th gradient term with respect to θ is an outer product:
Define wᵢ = σᵢ(1−σᵢ) and W = diag(w₁, …, wₙ). Summing the outer products over all i and writing it in matrix form:
Why this proves convexity. By section 2.2.1's Hessian characterization of convexity, L is convex everywhere its Hessian is positive semidefinite everywhere. Take any vector v:
Every term in that sum is a nonnegative weight wᵢ = σᵢ(1−σᵢ) ∈ (0, 1/4] times a squared number, so the whole sum is ≥ 0 for every v and every θ — H ⪰ 0 everywhere, with no exceptions and no dependence on where you are in parameter space. Binary cross-entropy is therefore convex on the whole of ℝ^{d+1}, which by 2.2.1 means: any point with zero gradient is automatically the global minimum, and there are no separate local minima anywhere else on the surface for gradient descent (or any other local method) to get trapped in.
When every wᵢ is strictly positive (true for any finite θ, since σᵢ is then strictly between 0 and 1) and X has full column rank, H is in fact strictly positive definite, giving a unique minimizer — not just a flat plateau of equally-good solutions. That strict positivity is exactly what breaks down as θ is driven toward infinity on separable data, which the Expert note below covers directly.
Where this is used: convexity guarantees a unique minimum exists, but unlike OLS — whose squared-error loss is also convex, and which gets a closed-form answer straight from the normal equations (section 2.3.1) with no iteration at all — logistic regression's loss involves the nonlinear σ function transcendentally, and setting ∇L(θ) = 0 has no algebraic solution for θ. Reaching that guaranteed minimum still requires an iterative method: plain gradient descent (as in the diagrams below), or a second-order method exploiting exactly this Hessian — Newton's method or IRLS, covered in sections 2.2.4 and 2.3.9.
Full-batch gradient descent on the cross-entropy loss, using the gradient ∇L = Xᵀ(σ(Xθ)−y) derived above. The boundary starts pointing the wrong way; scrub the slider or hit replay to watch it rotate into place as the loss falls.
The loss surface over two of the three parameters (bias fixed). Drag the green dot anywhere, release it, or pick a preset starting corner -- gradient descent always flows downhill to the same single minimum, exactly as the positive semidefinite Hessian above guarantees.
All three implementations fit the same model to the same kind of two-blob synthetic data: a bias plus two feature weights, minimizing binary cross-entropy. The first two hand-roll the exact gradient derived above; the third swaps in a production solver that reaches the same unregularized maximum-likelihood solution by a different (second-order) route.
#include <cmath>
#include <iostream>
#include <vector>
double sigmoid(double z) { return 1.0 / (1.0 + std::exp(-z)); }
double negLogLikelihood(const std::vector<double>& theta,
const std::vector<std::vector<double>>& X,
const std::vector<double>& y) {
double sum = 0.0;
for (size_t i = 0; i < X.size(); ++i) {
double z = theta[0] * X[i][0] + theta[1] * X[i][1] + theta[2] * X[i][2];
double p = sigmoid(z);
double eps = 1e-9;
sum += -(y[i] * std::log(p + eps) + (1 - y[i]) * std::log(1 - p + eps));
}
return sum / X.size();
}
std::vector<double> gradient(const std::vector<double>& theta,
const std::vector<std::vector<double>>& X,
const std::vector<double>& y) {
// The closed form derived above: Xᵀ(σ(Xθ) − y), averaged over n.
std::vector<double> g(3, 0.0);
for (size_t i = 0; i < X.size(); ++i) {
double z = theta[0] * X[i][0] + theta[1] * X[i][1] + theta[2] * X[i][2];
double err = sigmoid(z) - y[i];
g[0] += X[i][0] * err;
g[1] += X[i][1] * err;
g[2] += X[i][2] * err;
}
for (double& v : g) v /= X.size();
return g;
}
int main() {
// Small hand-built two-blob dataset with a bias column.
std::vector<std::vector<double>> X = {
{1, -2.0, -1.1}, {1, -1.1, -0.6}, {1, -1.6, -1.4}, {1, -0.8, -0.5},
{1, 1.2, 0.9}, {1, 1.8, 1.3}, {1, 0.9, 1.6}, {1, 2.1, 0.8},
};
std::vector<double> y = {0, 0, 0, 0, 1, 1, 1, 1};
std::vector<double> theta = {0.0, 0.0, 0.0};
double lr = 0.5;
for (int it = 0; it < 300; ++it) {
std::vector<double> g = gradient(theta, X, y);
for (int j = 0; j < 3; ++j) theta[j] -= lr * g[j];
if (it % 50 == 0) {
std::cout << "iter " << it << " loss=" << negLogLikelihood(theta, X, y) << "\n";
}
}
std::cout << "fitted theta: " << theta[0] << " " << theta[1] << " " << theta[2] << "\n";
return 0;
}- Spam detection. One of the original and still-common applications — word/token features feeding a linear score, squashed into "probability this is spam."
- Medical diagnosis screening. Predicting probability of a condition from clinical measurements, where the calibrated probability output (not just a hard yes/no) matters for setting a clinically appropriate decision threshold.
- Credit default / underwriting. Estimating the probability a borrower defaults from financial features — valued in regulated settings partly because coefficients are directly interpretable as log-odds effects (see the Pitfall below).
- Click-through-rate baselines. Before deep models, and still as a fast, well-understood baseline layer within them, logistic regression on engineered features is a standard first model for "will this ad/result be clicked."
- Churn prediction. Estimating the probability a customer cancels next period from usage and engagement features, ranking accounts by risk for retention outreach.
- Using squared error instead of cross-entropy "because it worked for OLS." As the Derivation above shows, squared error on a sigmoid output isn't the negative log-likelihood of any probabilistic model here, and its gradient carries a saturating
σ'factor that vanishes on confidently-wrong points — exactly the points training needs to correct fastest. - Reading a fitted coefficient
θⱼdirectly as "the change in probability per unit ofxⱼ." It isn't — it's the change in log-odds. Becauseσis nonlinear, the same coefficient produces a large probability swing nearp = 0.5and a tiny one nearp = 0orp = 1; onlye^θⱼas a multiplicative effect on the odds is coefficient-independent of where you evaluate it. - Hitting complete or quasi-complete separation — when some direction in feature space perfectly (or almost perfectly) separates the classes — and being surprised the fitted coefficients are enormous or the optimizer never stops improving. This isn't a bug; it's the Expert note below, and it needs a structural fix (regularization, or stopping early), not a smaller learning rate.
- Skipping feature scaling before gradient descent. The Hessian's eigenvalue spread (hence gradient descent's conditioning) depends directly on the relative scale of the columns of
X, exactly as it does for OLS — an unscaled feature with a much larger range dominates the step direction and slows convergence on every other feature.
Going deeper
Why separable data breaks maximum likelihood. Suppose the training data is linearly separable: some θ exists with yᵢ(θᵀxᵢ) > 0 for every example — every point strictly on the correct side of its boundary. Now consider scaling that same θ by a growing factor t. Every zᵢ = tθᵀxᵢ grows in magnitude in the correct direction, every σ(zᵢ) gets pushed strictly closer to the correct label, and the cross-entropy loss keeps strictly decreasing as t → ∞ — approaching 0 but never reaching it at any finite θ. The maximum likelihood estimate doesn't exist as a finite point; the optimizer just keeps inflating ‖θ‖ forever, chasing a supremum it can never attain. This is precisely where the strict positive-definiteness argument at the end of the Hessian derivation above breaks down: as θ grows this way, every wᵢ = σᵢ(1-σᵢ) → 0, so H itself degenerates toward the zero matrix along that direction — the loss surface goes flat, not because the model stopped learning, but because there's nothing left to distinguish. In practice this shows up as gradient descent never converging (loss keeps inching down, coefficients keep growing) or — in packages that use Newton-type solvers like IRLS (section 2.3.9) — as an explicit numerical warning about fitted probabilities collapsing to exactly 0 or 1. The standard fixes are exactly the ones Module 3 already built: an L2 penalty (ridge, section 2.3.2) or an L1 penalty (lasso, section 2.3.3) added to the loss keeps ‖θ‖ finite by directly penalizing its growth, and early stopping of the optimizer achieves something similar implicitly, without ever changing the loss function itself.
The Hessian H = XᵀWX (with W = diag(σᵢ(1−σᵢ))) is positive semidefinite for every θ, which proves the loss is convex everywhere. Why doesn't that alone guarantee a unique finite minimizer always exists?
Positive semidefinite only rules out negative curvature -- it guarantees no local minima exist apart from the global one, so any zero-gradient point found is the best point. But it says nothing about whether such a point exists at a finite θ. A unique minimizer additionally needs H to be strictly positive definite, which requires every weight wᵢ = σᵢ(1−σᵢ) to be strictly positive (true for any finite θ, since σᵢ is then strictly between 0 and 1) and X to have full column rank (no exactly collinear features). On linearly separable data, driving ‖θ‖ toward infinity in the separating direction pushes every σᵢ toward 0 or 1, so every wᵢ → 0 and H degenerates -- the loss keeps decreasing forever without ever reaching a finite minimum. Convexity guarantees there's only ever one basin to worry about; it doesn't guarantee that basin has a floor.
Logistic regression's whole story is one clean chain: squash a linear score with the sigmoid, read the model as making the log-odds linear in x, fit it by maximizing a Bernoulli likelihood, and get a gradient and a positive-semidefinite Hessian that together guarantee a single well-behaved minimum. Section 2.4.2, Softmax / Multinomial Regression, takes this exact same chain — one linear score per class, one shared normalization step in place of the sigmoid, one cross-entropy loss — and generalizes it from two classes to K.