Logistic Regression — Sigmoid, Cross-Entropy, from Scratch

Why this session exists

A yes/no decision often starts with a continuous score: how suspicious is a transaction, or how likely is a user to click? A linear predictor can return any real number, while a probability must lie between zero and one. Logistic regression connects those spaces through a model of log-odds. It is a useful baseline because its training objective is convex, its prediction is inexpensive, and its coefficients have a precise conditional interpretation—not because it is automatically calibrated or always the best classifier.

This lesson develops the probability model, derives its loss and gradient, implements training, and then separates three questions: does the model rank examples well, are its probabilities calibrated, and which decisions should we take? Those questions remain distinct in neural classifiers and tree ensembles.

By the end, you should be able to invert the logit, derive mean binary cross-entropy from a Bernoulli likelihood, implement stable gradient descent, select a threshold on validation data, and interpret a reliability table without mistaking a few observations for a calibration guarantee.

Prerequisites

Use the train/validation/test distinction from the ML mental model, matrix gradients from linear regression, and Bernoulli probabilities. Here an augmented design matrix $X\in\mathbb R^{n\times(d+1)}$ has a leading column of ones; $\theta$ includes the intercept. Labels are exactly 0 or 1. A dot product produces one logit per row.

(a) Intuition · 5 min

Think of a dimmer rather than an on/off switch. The linear score turns the dimmer; a separate threshold turns brightness into a binary action. A score of zero maps to probability 0.5, not necessarily to “take no action.” If missed positives are costly, acting at a much lower probability can be rational.

The three ideas to own before we code

1. Model: $p_i=\sigma(x_i^T\theta)$, with $\sigma(z)=1/(1+e^{-z})$.
2. Fit: minimize the mean negative Bernoulli log-likelihood. This estimates coefficients; it does not choose a business policy.
3. Decide: compare a score with a threshold selected for the intended population, costs and constraints. The same fitted model can support different operating points.

The sigmoid is not the only valid binary probability link. Probit uses a normal cumulative distribution function. Logistic regression specifically assumes that log-odds are linear; conditional on that assumption, the sigmoid is forced by algebra. Saying every binary classifier must use sigmoid confuses a model choice with a theorem.

A brief history

The useful connection to later models is mathematical rather than a company origin story. A neural network with no hidden layer, an affine scalar output and a sigmoid implements this same binary model. Adding hidden layers learns nonlinear features before the output; multinomial logistic regression instead uses a vector of logits and softmax for mutually exclusive classes. The scikit-learn 1.7 linear-model guide documents binary and multinomial formulations. Neither connection implies that every classifier or recommender has this architecture.

The sigmoid — the shape you will see forever

Start with odds $p/(1-p)$ and assume

Exponentiate: $p/(1-p)=e^z$. Rearrange $p=e^z(1-p)$ to get $p=e^z/(1+e^z)=\sigma(z)$. The inverse maps finite logits to the open interval $(0,1)$ mathematically; floating-point arithmetic can round the extremes to 0 or 1.

For orientation, $\sigma(0)=0.5$, $\sigma(2)\approx0.881$ and $\sigma(-2)\approx0.119$. Differentiating gives $\sigma'(z)=p(1-p)$, maximal at $p=0.5$. A fixed logit change therefore does not produce a fixed probability change.

Holding other represented features fixed, increasing feature $j$ by one unit multiplies the odds by $e^{\theta_j}$. A coefficient 0.7 gives an odds ratio about 2.014. Starting at $p=0.2$ gives odds 0.25; multiplying by 2.014 gives probability about 0.335, not 0.9 or a constant 0.7-point increase. Correlated predictors, feature transformations and confounding complicate interpretation; a coefficient is not a causal effect.

At threshold $t\in(0,1)$, predict positive when $x^T\theta\ge\log(t/(1-t))$. This moves a hyperplane in the represented feature space. Subtracting that cutoff from the decision score is equivalent for labels to shifting the intercept, but actually changing the fitted intercept changes the reported probabilities. Keep policy and probability estimation separately versioned.

Cross-entropy — where it comes from

For independent observations conditional on their inputs, the Bernoulli likelihood and its mean negative logarithm are

The three essential moves are write the Bernoulli probability, multiply across observations, take minus the mean log. Independence is a modeling assumption; repeated patients or temporal dependence also affect how we evaluate the fitted model.

For one example,

Multiplication cancels the denominator, so $\partial\ell/\partial z=p-y$ and

The Hessian is positive semidefinite: for any $v$, $v^TX^TWXv=\sum_iW_{ii}(x_i^Tv)^2\ge0$. This proves convexity, not strict convexity, existence of a finite minimizer, or convergence for every step size. Perfect separation can send unregularized weights toward infinity while loss approaches its infimum. L2 on slopes helps control that behavior when both classes are present; all-one-class data and unpenalized parameters need separate care.

For mean half squared error with linear predictions, the gradient also has form $X^T(\hat y-y)/n$. Without the half, it has an extra factor two. For squared error on sigmoid probabilities, an extra $p(1-p)$ remains; the parameter objective need not be convex. Squared probability error is nevertheless the meaningful Brier score, a proper scoring rule—not a nonsensical loss. BCE is the likelihood choice used here.

For numerical work use the logit expression

A confidently wrong logit $z=1000,y=0$ has loss about 1000, rather than infinity. Clipping probabilities can be useful for reporting but caps the loss; it is not the same as a stable logits implementation. Also, np.where evaluates both array arguments before selecting: writing two exponential formulas inside it does not avoid overflow. Evaluate only the relevant masked branch.

Metrics beyond accuracy — the whole zoo

Accuracy is the fraction of correct decisions at one threshold. Compare it with the majority-class baseline: at 0.1% prevalence, always-negative accuracy is 99.9% but positive recall is zero. Accuracy still measures something real; alone it does not describe rare-positive detection.

Precision is $TP/(TP+FP)$; recall is $TP/(TP+FN)$. F1 is $2TP/(2TP+FP+FN)$ when defined. State the positive label and zero-denominator convention. In the bundled breast-cancer dataset, the original target 1 means benign. The lab deliberately recodes malignant as 1 so recall refers to that class; it is an educational dataset, not a validated medical decision system.

ROC-AUC measures pairwise ranking, counting score ties as half a win. Its random-ranking reference is 0.5. Average precision (AP) summarizes precision weighted by recall increments; it is not in general trapezoidal PR area. The PR no-skill reference depends on positive prevalence. Neither AUC nor AP chooses a threshold, certifies calibration, or has a universal “shippable” range. See the metrics guide and the next metrics lesson.

Decision threshold — the knob nobody tunes

Lowering a fixed threshold can only add predicted positives, so recall cannot decrease when positives exist. Precision can go either way: admitting a true positive can improve it. Treat “lower threshold means lower precision” as a common tendency, not a monotonicity theorem. A threshold 0.5 does not balance precision and recall.

For calibrated deployment probabilities, zero cost for correct decisions, false-positive cost $C_{FP}$ and false-negative cost $C_{FN}$, the positive action costs $C_{FP}(1-p)$ in expectation, while the negative action costs $C_{FN}p$. Act positive if

Equal error costs give 0.5 without requiring balanced classes; prevalence is already inside calibrated posterior probabilities. Costs 5 and 500 give $5/505\approx0.00990$, not a guessed 0.1. If every review costs 5, including true positives, and review prevents the full 500 loss, the comparison instead gives $500p\ge5$, or 0.01. Capacity limits and imperfect interventions change the policy again.

For an empirical F1 or recall target, select the threshold on validation data and freeze it before testing. If repeatedly changing the policy after seeing test results, that set has become another validation set.

Calibration — reliability is a separate measurement

Group scores into bins. For each bin report mean predicted probability, fraction actually positive and sample count. Plot the first two against one another; the diagonal represents calibration. A bin with mean score 0.8 and 80 positives in 100 cases is consistent with calibration but uncertain; eight positives in ten cases is much weaker evidence. Empty bins have no estimate and should not be shown as zero-frequency observations.

Logistic regression's log-loss objective and canonical link can support good calibration under suitable specification and regularization, but do not guarantee it. Misspecified features, heavy penalties, class weighting, case-control sampling and distribution shift can distort the scale. An increasing score transform preserves ranking while altering calibration. Use held-out or cross-validated predictions for calibration fitting, and independent evaluation afterward. The calibration guide distinguishes reliability diagrams, sigmoid calibration and isotonic calibration; isotonic is flexible but needs enough calibration data.

(c) Hands-on · 25 min

This complete CPU lab uses a bundled dataset—no downloads. It separates training, validation and final test data, checks the gradient by finite differences, compares with scikit-learn, and reports reliability data for a plot. It prints measured results rather than promising a particular accuracy. It uses the scikit-learn 1.7 API.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (accuracy_score, roc_auc_score,
    average_precision_score, precision_recall_curve, confusion_matrix)
 
 
def sigmoid(z):
    z = np.asarray(z, dtype=float)
    out = np.empty_like(z)
    pos = z >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-z[pos]))
    ez = np.exp(z[~pos])
    out[~pos] = ez / (1.0 + ez)
    return out
 
 
def objective(theta, X, y, lam=0.01):
    z = X @ theta
    penalty = theta.copy()
    penalty[0] = 0  # do not penalize the intercept
    loss = np.mean(np.logaddexp(0, z) - y*z)
    loss += 0.5 * lam * (penalty @ penalty)
    grad = X.T @ (sigmoid(z) - y) / len(y) + lam * penalty
    return loss, grad
 
 
def fit_logreg(X, y, lr=0.1, epochs=1800, lam=0.01):
    theta = np.zeros(X.shape[1])
    initial = objective(theta, X, y, lam)[0]
    for _ in range(epochs):
        loss, grad = objective(theta, X, y, lam)
        theta -= lr * grad
    final, grad = objective(theta, X, y, lam)
    assert np.isfinite(final) and final < initial
    print('loss', initial, '->', final, 'gradient norm', np.linalg.norm(grad))
    return theta
 
 
def reliability(y, p, bins=8):
    idx = np.minimum((p*bins).astype(int), bins-1)
    rows = []
    for b in range(bins):
        mask = idx == b
        if mask.any():
            rows.append((float(p[mask].mean()), float(y[mask].mean()), int(mask.sum())))
    return rows
 
 
X, original_y = load_breast_cancer(return_X_y=True)
y = 1 - original_y  # positive means malignant
Xdev, Xte, ydev, yte = train_test_split(X, y, test_size=.2, stratify=y, random_state=42)
Xtr, Xva, ytr, yva = train_test_split(Xdev, ydev, test_size=.25, stratify=ydev, random_state=43)
scaler = StandardScaler().fit(Xtr)
add_bias = lambda A: np.column_stack([np.ones(len(A)), scaler.transform(A)])
A, V, T = (add_bias(Z) for Z in (Xtr, Xva, Xte))
 
# Independent numerical oracle for every coordinate of the analytic gradient.
rng = np.random.default_rng(5)
test_theta = rng.normal(scale=.1, size=A.shape[1])
eps = 1e-6
analytic = objective(test_theta, A[:12], ytr[:12])[1]
numeric = np.empty_like(analytic)
for j in range(len(test_theta)):
    delta = np.zeros_like(test_theta); delta[j] = eps
    numeric[j] = (objective(test_theta+delta, A[:12], ytr[:12])[0]
                  - objective(test_theta-delta, A[:12], ytr[:12])[0])/(2*eps)
assert np.allclose(analytic, numeric, atol=1e-7)
with np.errstate(over='raise', invalid='raise', divide='raise'):
    assert np.allclose(sigmoid(np.array([-1000., 0., 1000.])), [0., .5, 1.])
 
theta = fit_logreg(A, ytr)
pva = sigmoid(V @ theta)
precision, recall, thresholds = precision_recall_curve(yva, pva)
f1 = np.divide(2*precision*recall, precision+recall,
               out=np.zeros_like(precision), where=(precision+recall)>0)
best = int(np.argmax(f1[:-1]))  # terminal PR point has no threshold
cutoff = float(thresholds[best])
print('validation cutoff', cutoff, 'validation F1', f1[best])
 
# Match sklearn's sum-loss + ||w||^2/(2C) convention to mean-loss + lam||w||^2/2.
sk = LogisticRegression(C=1/(len(ytr)*.01), max_iter=2000, tol=1e-9)
sk.fit(scaler.transform(Xtr), ytr)
print('validation probability difference', np.max(np.abs(pva-sk.predict_proba(scaler.transform(Xva))[:, 1])))
 
# Final test: no more choices made from these results.
pte = sigmoid(T @ theta)
for t in (.5, cutoff):
    print('threshold', t, 'accuracy', accuracy_score(yte, pte >= t))
    print(confusion_matrix(yte, pte >= t, labels=[0, 1]))
print('test ROC-AUC', roc_auc_score(yte, pte), 'AP', average_precision_score(yte, pte))
rows = reliability(yte, pte)
assert sum(r[2] for r in rows) == len(yte)
print('reliability rows (mean score, positive fraction, count):', rows)

What each block does

The scaler learns only training means and variances. add_bias adds the intercept after scaling. objective uses the same mean-loss convention for its scalar and gradient, with a separately stated L2 penalty; the intercept is exempt. Gradient differences test the implementation locally, not its convergence on every dataset. The optimizer reports its final gradient norm so a finite epoch budget is not silently called an exact optimum.

The scikit-learn comparison matches the regularization convention as well as the data. L-BFGS is a quasi-Newton solver; finite-step gradient descent can differ from it even with the same objective. Compare loss and predictions, not a promised number of matching decimal places.

precision_recall_curve uses positive predictions at scores greater than or equal to its thresholds; the lab uses the same comparison. Its last precision/recall pair has no corresponding threshold and is excluded when selecting one. The test reliability table is a diagnostic, not a calibration fit. To draw it, put each row's mean score on x and positive fraction on y, annotate counts, and add the diagonal.

(d) Production reality · 15 min

Linear model, engineered features, or trees?

A plain logistic model is additive in log-odds. It is attractive when that assumption is useful, when a compact model matters, or when you need an inspectable baseline. Convexity removes bad local minima in the stated objective, but not collinearity, data leakage, unstable coefficients or an unsuitable representation. Standardization also changes what an equal penalty on coefficients means.

For an interaction such as “high spending is suspicious only for newly created accounts,” add a spending-by-account-age feature or a suitable basis expansion. Splines and interaction terms preserve convex optimization in the coefficients, but their choice is substantive modeling work. Feature engineering can make a logistic decision boundary nonlinear in the raw inputs.

Boosted trees can learn threshold interactions automatically and may improve validation performance. They also bring tree-count/depth choices, a larger serving artifact and their own calibration checks. Missing-value and categorical support depend on the implementation. Neither model family is automatically compliant, causal or clinically safe; interpretability is one property to assess with the actual review requirements. Keep a model only when its measured benefit justifies the complexity, rather than claiming a company universally uses it.

Combining an explicit linear branch with learned features

A hybrid can add a linear score on engineered features to a learned nonlinear score before the sigmoid:

Here $\phi(x)$ might contain explicit category crosses that memorize frequently observed combinations, while $g_\psi$ learns a representation shared across inputs. The final BCE derivative with respect to the combined logit remains $p-y$, but fitting $\psi$ makes the overall parameter problem nonconvex. This is a mathematical example of combining linear and neural mechanisms, not a claim that any particular recommender architecture is identical to this formula. Both branches need the same feature-availability boundary, and adding the nonlinear branch does not automatically improve ranking or calibration.

Where this shows up in the rest of the plan

Regularization explains the L2 term and alternatives. Bias and variance distinguishes approximation error from variation across fits. Decision trees provide a nonlinear comparison. Evaluation metrics develops the ROC/AP and capacity distinctions. The same scalar sigmoid can appear after a neural feature extractor, but softmax, multilabel sigmoid heads and other objectives require their own contracts.

(e) Recall + stretch · 10 min

1. Invert the logit. Answer: exponentiate $\log(p/(1-p))=z$, then solve $p=e^z(1-p)$ to obtain sigmoid. The assumption is linear log-odds, not merely “binary labels.”
2. Write mean BCE and its gradient. Answer: the negative mean Bernoulli log-likelihood above; $X^T(p-y)/n$. The L2 lab adds $\lambda\theta$ only on slopes.
3. Why not blindly substitute MSE? Answer: squared probability error is meaningful, but with a sigmoid it changes the objective and adds the potentially saturating derivative. It is not the Bernoulli maximum-likelihood calculation.
4. What does 99% accuracy conceal? Answer: at 1% prevalence the always-negative classifier already gets 99%. Its positive recall is 0, and its precision is undefined unless a reporting convention is supplied.
5. ROC-AUC or AP? Answer: ROC-AUC describes pairwise ranking; AP weights precisions at recall increments and depends on prevalence. Report an operating-point confusion matrix too.
6. Where is the threshold selected? Answer: validation predictions, not final test labels. Select for a declared objective, freeze, then evaluate uncertainty on test.
7. Does probability 0.8 mean 80% accuracy? Answer: calibration asks whether similarly scored examples have about 80% positive outcomes; classification accuracy uses thresholded labels and is different.
8. Why might a linear model be preferred to a booster? Answer: simpler serving, interpretable conditional coefficients, adequate held-out quality or maintenance cost. No legal or universal performance conclusion follows from the model name.
9. Imbalance experiment. Downsample malignant training cases only, keeping validation/test representative. Compare unweighted and class_weight='balanced' fits using malignant-positive recall, AP and reliability. Answer: weighting changes the fitted objective and often the probability scale; neither AP recovery nor a particular accuracy is guaranteed. Do not downsample the test set to manufacture a nicer precision.
10. Cost stretch. With false-positive cost 5 and false-negative cost 500, zero costs otherwise, derive $t=5/505$. Enumerate validation thresholds to compare empirical costs. If only 20 reviews are possible, also evaluate top-20 precision/recall; the unconstrained Bayes threshold may exceed capacity.

Explain-out-loud test

Explain the chain linear log-odds → sigmoid → Bernoulli likelihood → BCE → gradient, then explain why calibration and threshold selection are not consequences of that chain. A complete answer includes the log-odds assumption, the mean-loss normalization, and the data split used for each decision.

What comes next

Continue with Regularization — L1, L2, Elastic Net, where coefficient constraints become a deliberate bias–variance choice rather than an unexplained library default.

Sources and execution boundary

Reviewed against the scikit-learn 1.7 logistic formulation, calibration guide, PR API and average precision definition, accessed 2026-10-07. The bounded local CPU run checks this program and its numerical assertions. It is not a clinical validation, production benchmark, exhaustive optimization proof, or evidence about a particular company's deployment.