03 - Logistic Regression and Binary Classification

A hospital wants to estimate, from thirty measurements of a nodule, whether it is benign or malignant. The answer it needs is not "number 0.3" as a straight line would produce, but a probability between 0 and 1, to which somebody will then apply a threshold chosen according to the cost of the errors. Logistic regression does exactly this: it takes the linear combination of the measurements from chapter 02 and squeezes it into the interval [0, 1] with an S-shaped function, the sigmoid. Despite the name it is a classifier. Here we derive it, train it in two ways that reach the same minimum, compare it with scikit-learn, and see how to read its probabilities when the classes are imbalanced or more than two.

The examples use scikit-learn 1.9.0, numpy 2.4.6 and scipy 1.17.1, with the Breast Cancer Wisconsin (569 nodules, 30 measurements) and Digits (1797 8x8 images of digits) datasets bundled in the library, plus synthetic data with a fixed seed. The outputs are identical across two runs.

From the line to the probability

Let z = w·x + b be the linear score of an example. The sigmoid σ(z) = 1/(1 + e^−z) turns any real number into a value between 0 and 1: σ(0) = 0.5, very positive values give probabilities close to 1, very negative ones close to 0. The model says P(y = 1 | x) = σ(z). The meaning of the weights is best read in terms of odds, the ratio p/(1 − p) between the probability of an event and that of its opposite: a 75% probability is equivalent to odds of 3 to 1. Inverting the sigmoid gives z = ln(p/(1 − p)), the logarithm of the odds: the model is linear in the log-odds, and increasing a measurement by one unit multiplies the odds by e^w. The textbook by James, Witten, Hastie and Tibshirani (An Introduction to Statistical Learning, section 4.3) develops the same reading.

The decision boundary, that is, the place where the model is undecided (p = 0.5), is w·x + b = 0: a hyperplane. Logistic regression separates the classes with flat surfaces; if the true boundary is curved you need transformed measurements or another model, such as those in chapters 06 and 07.

The loss: cross-entropy

If each example is a coin with probability p_i = σ(z_i), the likelihood of the set is the product of the probabilities assigned to the observed labels, and it is convenient to minimise the opposite of its logarithm: L = −(1/n) Σ [y_i ln p_i + (1 − y_i) ln(1 − p_i)], called cross-entropy or log loss. The gradient with respect to the weights has a very simple form, (1/n) Xᵀ(p − y), the same as in linear regression with p in place of the prediction. The Hessian is (1/n) XᵀWX with W = diag(p_i(1 − p_i)), positive semidefinite: the loss is convex, so it has a single minimum and any method that descends finds it.

Worked exercise: one prediction, one loss

An example has score z = 0.8. The odds are e^0.8 = 2.2255 and the probability σ(0.8) = 0.69. If the true label is 1 the loss is −ln 0.69 = 0.3711; if it is 0 the loss is −ln 0.31 = 1.1711. The same prediction costs more than three times as much when it is wrong. The gradient with respect to z is p − y: −0.31 if y = 1 (it pushes z up) and +0.69 if y = 0 (it pushes it down, harder, because the error is larger).

import numpy as np
from scipy.special import expit

z = 0.8                                   # linear score w.x + b of one example
p = expit(z)
print("odds = e^z:", round(float(np.exp(z)), 4), " p = sigmoid(z):", round(float(p), 4))
print("loss if the true label is 1:", round(float(-np.log(p)), 4))
print("loss if the true label is 0:", round(float(-np.log(1 - p)), 4))
print("gradient with respect to z (p - y), y=1 and y=0:", round(float(p - 1), 4), round(float(p - 0), 4))

# Output:
# odds = e^z: 2.2255  p = sigmoid(z): 0.69
# loss if the true label is 1: 0.3711
# loss if the true label is 0: 1.1711
# gradient with respect to z (p - y), y=1 and y=0: -0.31 0.69

The printed values confirm the hand computations. The loss grows without limit when the model is sure and wrong (−ln p tends to infinity for p → 0): this is why cross-entropy punishes wrong confidence far more than a moderate error, and why in programs it is computed with np.logaddexp(0, z) − y·z instead of with log(sigmoid), which for very negative z would produce log(0).

Newton-Raphson and gradient descent

The listing implements the loss with an L2 penalty, its gradient and its Hessian, checks the gradient with finite differences, and minimises in two ways. Gradient descent uses a fixed step of 0.5. The Newton-Raphson method (called IRLS in statistics, iteratively reweighted least squares) at each step solves H·δ = g and moves the weights by δ: it uses the curvature to choose the length of the step, and near the minimum it converges in very few iterations. The penalty is written in the scikit-learn convention: in the documentation, the minimised function is the mean of the loss plus r(W)/(S·C), with r = ½||W||² and S equal to the number of examples when there are no sample weights. With C = 1 this is equivalent to my lam = 1/(C·n), applied to the weights and not to the intercept.

import numpy as np
from scipy.special import expit
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)
sc = StandardScaler().fit(Xtr)
A = np.c_[np.ones(len(Xtr)), sc.transform(Xtr)]
n, d = A.shape
C = 1.0
lam = 1 / (C * n)                    # same penalty as sklearn: mean loss + (lam/2)||w||^2, intercept excluded

def loss(w):
    z = A @ w
    ll = np.mean(np.logaddexp(0, z) - ytr * z)           # cross-entropy, stable for large z
    return ll + 0.5 * lam * np.sum(w[1:] ** 2)

def grad(w):
    g = A.T @ (expit(A @ w) - ytr) / n
    g[1:] += lam * w[1:]
    return g

def hess(w):
    p = expit(A @ w)
    H = (A * (p * (1 - p))[:, None]).T @ A / n
    H[1:, 1:] += lam * np.eye(d - 1)
    return H

# 1) gradient check with finite differences
rng = np.random.default_rng(0)
w = 0.1 * rng.standard_normal(d)
num = np.array([(loss(w + 1e-6 * e) - loss(w - 1e-6 * e)) / 2e-6 for e in np.eye(d)])
print("analytic gradient = finite differences:", np.allclose(num, grad(w), atol=1e-8))

# 2) Newton-Raphson (IRLS)
w = np.zeros(d)
for it in range(1, 30):
    step = np.linalg.solve(hess(w), grad(w))
    w = w - step
    if np.linalg.norm(step) < 1e-10:
        break
print("Newton: iterations =", it, " loss =", round(loss(w), 6))

# 3) gradient descent with a fixed step
wg = np.zeros(d)
for k in range(1, 20001):
    wg = wg - 0.5 * grad(wg)
    if np.linalg.norm(grad(wg)) < 1e-6:
        break
print("gradient descent: iterations =", k, " loss =", round(loss(wg), 6))

# 4) comparison with scikit-learn (same C, same L2 penalty)
sk = LogisticRegression(C=C, max_iter=5000, tol=1e-10).fit(sc.transform(Xtr), ytr)
w_sk = np.r_[sk.intercept_, sk.coef_.ravel()]
print("Newton = sklearn (tolerance 1e-5):", np.allclose(w, w_sk, atol=1e-5))
print("gradient descent = Newton (tolerance 1e-3):", np.allclose(wg, w, atol=1e-3))
At = np.c_[np.ones(len(Xte)), sc.transform(Xte)]
print("test accuracy:", round(float(np.mean((expit(At @ w) > 0.5) == yte)), 4))

# 5) without a penalty (C=inf) the whole dataset is linearly separable: the weights have no maximum
Xs = StandardScaler().fit_transform(X)
for c in (1.0, np.inf):
    m = LogisticRegression(C=c, max_iter=5000, tol=1e-10).fit(Xs, y)
    print(f"C={c}: largest absolute weight = {np.abs(m.coef_).max():.1f}, train accuracy = {m.score(Xs, y):.3f}")

# Output:
# analytic gradient = finite differences: True
# Newton: iterations = 10  loss = 0.068116
# gradient descent: iterations = 5963  loss = 0.068116
# Newton = sklearn (tolerance 1e-5): True
# gradient descent = Newton (tolerance 1e-3): True
# test accuracy: 0.986
# C=1.0: largest absolute weight = 1.3, train accuracy = 0.988
# C=inf: largest absolute weight = 64176.2, train accuracy = 1.000

The analytic gradient coincides with the numerical one. Newton reaches the minimum (loss 0.068116) in 10 iterations; gradient descent takes 5963 to reach the same loss, because it does not know the curvature. Both coincide with LogisticRegression(C=1.0) within the stated tolerances. Newton costs more at each step, because it solves a d × d system, and it is not convenient when the measurements number in the thousands; for this scikit-learn offers lbfgs (the default), which estimates the curvature from past gradients.

The test accuracy is 0.986, that is, 141 correct cases out of 143. Removing the penalty (C=np.inf; the penalty argument is deprecated in the installed version and C and l1_ratio are used instead) changes the behaviour: in the last lines of the listing, on the whole dataset, the data turn out to be separable by a hyperplane, the largest weight exceeds 60000 (a value that depends on the solver tolerance: it has no limit) and the training accuracy is 1.000, while with C=1.0 the largest weight stays below 1.5. When the classes are separable the loss has no finite minimum, because it is enough to let the weights grow to make the probabilities ever sharper; the penalty serves this purpose too.

From probability to decision

The model returns probabilities; the decision requires a threshold. The terms are those of the confusion matrix: true positives (TP), false positives (FP), false negatives (FN), true negatives (TN). Precision is TP/(TP+FP), the share of correct alarms; recall is TP/(TP+FN), the share of positives found. ROC-AUC and PR-AUC summarise the quality of the ranking over all thresholds, and the log loss and the Brier score (mean squared error of the probabilities) measure the probabilities themselves. On the dataset the positive class is now the malignant tumour (53 cases out of 143 in the test set).

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (average_precision_score, brier_score_loss, confusion_matrix,
                             log_loss, precision_score, recall_score, roc_auc_score)
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
y = 1 - y                                   # in the dataset 0 = malignant: we make it the "positive" class
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)).fit(Xtr, ytr)
p = model.predict_proba(Xte)[:, 1]
print("ROC-AUC:", round(roc_auc_score(yte, p), 4), " PR-AUC:", round(average_precision_score(yte, p), 4))
print("log loss:", round(log_loss(yte, p), 4), " Brier:", round(brier_score_loss(yte, p), 4))
for threshold in (0.5, 0.1):
    pred = (p >= threshold).astype(int)
    tn, fp, fn, tp = confusion_matrix(yte, pred).ravel()
    print(f"threshold {threshold}: TP={tp} FP={fp} FN={fn} TN={tn}  precision={precision_score(yte, pred):.3f}  recall={recall_score(yte, pred):.3f}")

# Output:
# ROC-AUC: 0.9962  PR-AUC: 0.9943
# log loss: 0.0776  Brier: 0.0224
# threshold 0.5: TP=49 FP=1 FN=4 TN=89  precision=0.980  recall=0.925
# threshold 0.1: TP=52 FP=6 FN=1 TN=84  precision=0.897  recall=0.981

With threshold 0.5 the model finds 49 malignant cases out of 53 (recall 0.925) with a single false alarm; lowering the threshold to 0.1 the missed malignant cases drop to 1 (recall 0.981) and the false alarms rise to 6, the precision falls to 0.897. Neither threshold is "right": it depends on the relative cost of missing a tumour and of sending a person to one more examination. ROC-AUC 0.9962 and PR-AUC 0.9943 say that the ranking is excellent, not which threshold to use. With 53 positives in the test set, a single case moves the recall by almost two points: these figures are to be read as orders of magnitude and not as a ranking between models. The full family of metrics is in chapter 11.

Imbalanced classes

When 97% of the examples are negative, a model that always says "negative" has 97% accuracy. We build a synthetic case with 3% positives and look at three choices: no intervention, lowering the threshold, and the parameter class_weight="balanced", which multiplies the loss of each example by a weight inversely proportional to the frequency of its class.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import average_precision_score, precision_score, recall_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=5000, n_features=10, n_informative=4, weights=[0.97, 0.03],
                           class_sep=1.5, flip_y=0.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
print("share of positives in the test set:", round(yte.mean(), 3), " accuracy of 'always 0':", round(1 - yte.mean(), 3))

def report(name, pred, p):
    print(f"{name:34s} precision={precision_score(yte, pred, zero_division=0):.3f} "
          f"recall={recall_score(yte, pred):.3f} PR-AUC={average_precision_score(yte, p):.3f}")

plain = LogisticRegression(max_iter=1000).fit(Xtr, ytr)
p0 = plain.predict_proba(Xte)[:, 1]
report("no weights, threshold 0.5", (p0 >= 0.5).astype(int), p0)
report("no weights, threshold 0.15", (p0 >= 0.15).astype(int), p0)
bal = LogisticRegression(class_weight="balanced", max_iter=1000).fit(Xtr, ytr)
p1 = bal.predict_proba(Xte)[:, 1]
report("class_weight='balanced', threshold 0.5", (p1 >= 0.5).astype(int), p1)
print("mean predicted probability: no weights", round(p0.mean(), 3), " balanced", round(p1.mean(), 3))

# Output:
# share of positives in the test set: 0.03  accuracy of 'always 0': 0.97
# no weights, threshold 0.5          precision=0.583 recall=0.311 PR-AUC=0.354
# no weights, threshold 0.15         precision=0.383 recall=0.511 PR-AUC=0.354
# class_weight='balanced', threshold 0.5 precision=0.110 recall=0.733 PR-AUC=0.170
# mean predicted probability: no weights 0.035  balanced 0.295

Without weights and with threshold 0.5 the model is cautious: precision 0.583 and recall 0.311. Lowering the threshold to 0.15 shifts the balance (recall 0.511, precision 0.383) without changing the model, and leaves the PR-AUC at 0.354. With class_weight="balanced" the recall rises to 0.733 but the precision collapses to 0.110, and the PR-AUC is 0.170, that is, the ranking got worse; moreover the mean predicted probability goes from 0.035 (close to the true frequency, 0.03) to 0.295, so the probabilities are no longer probabilities.

The result of this experiment, which holds for these synthetic data and not for every problem, is that weighting the classes changes the model and not only the threshold; if all you need is more recall, moving the threshold is cheaper and leaves the probabilities interpretable. The scikit-learn documentation on calibration observes that logistic regression usually returns well-calibrated probabilities because it optimises the log loss, and class weights break exactly that property.

More than two classes: softmax

With K classes there are K weight vectors and the probability of class k is the softmax of the scores: exp(z_k) / Σ_l exp(z_l), positive and summing to 1 by construction. The alternative is one-vs-rest: one binary classifier per class, which does not guarantee probabilities that sum to 1. In version 1.9.0 the multi_class argument of LogisticRegression no longer exists (passing it raises a TypeError): the multiclass case is multinomial by default, and one-vs-rest is obtained by wrapping the model in OneVsRestClassifier.

import numpy as np
from scipy.special import softmax
from sklearn.datasets import load_digits
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.multiclass import OneVsRestClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_digits(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=0, stratify=y)
soft = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(Xtr, ytr)
ovr = make_pipeline(StandardScaler(), OneVsRestClassifier(LogisticRegression(max_iter=2000))).fit(Xtr, ytr)
print("8x8 digits: classes, features:", len(set(y)), X.shape[1])
print("softmax (multinomial):", round(soft.score(Xte, yte), 4), " one-vs-rest:", round(ovr.score(Xte, yte), 4))

lr = soft[-1]
Z = lr.decision_function(soft[0].transform(Xte[:3]))
print("hand-made softmax = predict_proba:", np.allclose(softmax(Z, axis=1), soft.predict_proba(Xte[:3])))
print("sum per row:", softmax(Z, axis=1).sum(axis=1).round(6).tolist())

# Output:
# 8x8 digits: classes, features: 10 64
# softmax (multinomial): 0.9689  one-vs-rest: 0.9689
# hand-made softmax = predict_proba: True
# sum per row: [1.0, 1.0, 1.0]

On the 8x8 digits with 10 classes the two schemes have the same accuracy (0.9689): on this problem and with this seed there is no difference, and it cannot be attributed to any general rule. The softmax computed by hand from decision_function coincides with predict_proba, and the rows sum to 1.

When not to use logistic regression

• The boundary between the classes is not flat and you cannot build measurements that make it so.
• You have many interactions between measurements that you do not know in advance.
• The target is a number and not a category: you need regression.

In the first two cases a tree or a set of trees captures the structure without you specifying it, at the cost of explainability; in the third the problem is badly formulated. Where, instead, explaining every prediction matters, logistic regression has the advantage that each weight is an odds multiplier that can be read and discussed.

What you can do now

You can derive the loss and the gradient of logistic regression and check them, train it with Newton and with gradient descent and recognise the scikit-learn convention, choose a threshold after reading the confusion matrix, distrust class weights when you need probabilities, and use the softmax for more classes. The next chapter explains why a model can be too simple or too flexible: 04 - Bias-Variance Tradeoff.

The same series also includes: 10 - Feedforward Neural Networks: From Perceptron to Deep Learning and 12 - ML Deployment in Production: Serving and Monitoring.