Logistic Regression: A Complete Guide to Binary Classification

Logistic regression was one of the first models that helped me connect a machine-learning equation to a complete training process. The model is compact enough to calculate by hand, but it contains the same broad ingredients used by neural networks: a linear transformation, an activation function, a loss, gradients, and parameter updates.

This guide develops logistic regression for binary classification from the prediction equation through a complete NumPy implementation. Along the way, I explain what the parameters mean, why binary cross-entropy is used, how the gradients are derived, and how vectorization trains on many examples without an explicit loop over the dataset.

The binary-classification problem

Binary classification assigns an input to one of two classes. The target is commonly encoded as:

\[y\in\{0,1\}\]

For a cat-image classifier, \(y=1\) might mean that an image contains a cat, while \(y=0\) means that it does not. A \(64\times64\) RGB image contains three color values at each pixel, so flattening it creates:

\[n_x=64\times64\times3=12{,}288\]

The resulting input vector is \(x\in\mathbb{R}^{n_x}\). Logistic regression uses this vector to estimate the conditional probability of the positive class. The subscript below identifies the fitted model distribution; its estimate can differ from the true conditional probability:

\[\hat{y}=P_{w,b}(y=1\mid x)\]

For example, the same setup applies to spam detection, fraud screening, medical-test classification, and many other two-class problems. In every case, the meaning of class 1 must be defined consistently before training.

From a linear score to a probability

The model first calculates an unrestricted linear score:

\[z=w^Tx+b\]

• \(x\in\mathbb{R}^{n_x\times1}\) is the input.
• \(w\in\mathbb{R}^{n_x\times1}\) contains one weight per feature.
• \(b\in\mathbb{R}\) is the bias or intercept.
• \(w^Tx\) and \(z\) are scalars for one example.

The score can be any real number. Logistic regression maps it to a probability using the sigmoid function:

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

The model probability is:

\[\boxed{\hat{y}=\sigma(w^Tx+b)}\]

Sigmoid maps every finite real input to a value strictly between zero and one. When \(z=0\), the output is 0.5. Large positive scores approach one, and large negative scores approach zero:

\[z\to+\infty\Rightarrow\sigma(z)\to1,\qquad z\to-\infty\Rightarrow\sigma(z)\to0\]

For example, a prediction of 0.90 means that the fitted model assigns an estimated 90% probability to class 1. Calibration asks whether, among many cases assigned probabilities near 0.90, about 90% are positive. That must be checked on representative held-out data.

This linear-transformation-plus-activation pattern also appears in neural networks. See Neural Networks for Beginners: From Architecture to Forward Propagation for the extension from one logistic unit to multiple layers.

Weights, bias, and log-odds

Expanding the score makes the role of the parameters visible:

\[z=w_1x_1+w_2x_2+\cdots+w_{n_x}x_{n_x}+b\]

Holding other features constant, a positive \(w_j\) makes the positive class more likely as \(x_j\) increases; a negative value makes it less likely. More precisely, logistic regression is linear in the log-odds:

\[\log\left(\frac{\hat{y}}{1-\hat{y}}\right)=w^Tx+b\]

Odds are the ratio of the positive-class probability to the negative-class probability: a probability of 0.8 gives odds of \(0.8/0.2=4\), or 4:1. Increasing \(x_j\) by one unit changes the log-odds by \(w_j\), holding other variables constant, and multiplies the odds by \(e^{w_j}\). The change in probability depends on its starting value because sigmoid is nonlinear. These coefficients describe associations in the fitted model; they do not by themselves establish the effect of intervening on a feature.

The bias shifts the score. At threshold \(t=0.5\), a nonzero weight vector with no bias gives a decision boundary through the origin. The intercept allows that boundary to shift:

\[w^Tx+b=0\]

Classification threshold and decision boundary

A probability estimate can be converted into a label with a threshold \(t\):

\[\hat{y}_{\text{class}}=\begin{cases}1,&\hat{y}\ge t\\0,&\hat{y}<t\end{cases}\]

At the common default \(t=0.5\), sigmoid’s monotonicity makes this equivalent to predicting class 1 whenever \(w^Tx+b\ge0\). However, 0.5 is not mandatory. Fraud screening may favor recall and use a lower threshold, while a high-cost intervention may require a higher threshold. For any threshold \(0<t<1\), the boundary is \(w^Tx+b=\log(t/(1-t))\). Select the threshold on validation data and fix it before evaluating the final decision rule on test data.

With two input features and nonzero \(w\), the boundary is a straight line; in higher dimensions it is a hyperplane. Sigmoid changes the probability along the score axis but does not bend this boundary. Curved boundaries require transformed features or a more flexible model. For instance, adding \(x_1x_2\) as a feature lets logistic regression use interactions between the two original inputs.

Binary cross-entropy loss

Training needs an objective for choosing the parameters. The usual choice for logistic regression is binary cross-entropy, also called log loss; all logarithms here are natural:

\[\mathcal{L}(\hat{y},y)=-\left[y\log(\hat{y})+(1-y)\log(1-\hat{y})\right]\]

For a positive example, \(y=1\), the second term disappears:

\[\mathcal{L}(\hat{y},1)=-\log(\hat{y})\]

The loss is small when \(\hat{y}\) is close to one and becomes extremely large when the model confidently predicts a value near zero. In contrast, for a negative example, \(y=0\):

\[\mathcal{L}(\hat{y},0)=-\log(1-\hat{y})\]

Now the loss is small near zero and large near one. Cross-entropy therefore penalizes confident mistakes much more strongly than uncertain ones.

From individual loss to dataset cost

For \(m\) examples, the cost is the average loss:

\[J(w,b)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}\left(\hat{y}^{(i)},y^{(i)}\right)\]

Training searches for \(w\) and \(b\) that minimize this cost. A Bernoulli label has model probability \(a^y(1-a)^{1-y}\): this selects \(a\) when \(y=1\) and \(1-a\) when \(y=0\). Modeling the labels as independent conditional on their inputs makes the dataset likelihood the product of those terms. Taking its negative logarithm gives the sum of BCE losses, so minimizing their average gives the same maximum-likelihood fit.

Why not use squared error?

It is possible to combine sigmoid with squared error, but it is usually not the preferred formulation for logistic regression. With sigmoid, squared error can produce a less convenient, non-convex optimization objective and gradients that shrink when sigmoid saturates. Binary cross-entropy matches the Bernoulli probability model and yields a convex objective for ordinary logistic regression.

Squared error on a binary probability estimate is also a valid scoring rule: its conditional expected value is minimized by the true positive-class probability. The optimization difficulty above comes from composing that loss with sigmoid and fitting the parameters. For BCE, \(d^2\mathcal L/dz^2=a(1-a)\ge0\); composing this convex logit loss with the affine score and averaging preserves convexity. Convexity does not guarantee a unique or finite minimizer, as the separable dataset below illustrates.

Deriving the gradients

For one example, define \(a=\hat{y}=\sigma(z)\). The derivative of sigmoid is:

\[\frac{da}{dz}=a(1-a)\]

Differentiate the loss with respect to the probability, then multiply by the sigmoid derivative:

\[\frac{d\mathcal{L}}{dz}=\frac{d\mathcal{L}}{da}\frac{da}{dz}=\left(-\frac{y}{a}+\frac{1-y}{1-a}\right)a(1-a)=a-y\]

For example, \(a=0.5\) and \(y=1\) give \(d\mathcal{L}/dz=-0.5\). For a small calculation, hold \(w\) fixed and update only \(b\) with learning rate 0.1: \(b\leftarrow b-0.1(-0.5)=b+0.05\). This raises \(z\) from 0 to 0.05 and the probability from 0.5 to about 0.5125. Because \(z=w^Tx+b\), the parameter derivatives are:

\[\frac{\partial\mathcal{L}}{\partial w}=x(a-y),\qquad \frac{\partial\mathcal{L}}{\partial b}=a-y\]

Finally, averaging over all examples gives the cost gradients:

\[dw=\frac{1}{m}X(A-Y)^T\]

\[db=\frac{1}{m}\sum_{i=1}^{m}\left(a^{(i)}-y^{(i)}\right)\]

Here \(X\) collects the input columns, \(A\) collects their predicted probabilities, and \(Y\) collects their labels. Both \(A\) and \(Y\) have shape \((1,m)\). Thus \(X(A-Y)^T\) has shape \((n_x,1)\), matching \(w\). In the code below, dZ = A - Y stores the individual examples’ logit derivatives; the division by \(m\) occurs when forming dw and db. Using the gradient of the mean cost directly at dZ instead would require omitting those later divisions.

For a more general computational-graph treatment, continue with Derivatives and Computation Graphs for Neural Network Learning.

Gradient descent

Gradient descent moves the parameters in the direction opposite the gradient:

\[w:=w-\alpha dw\]

\[b:=b-\alpha db\]

Here, the learning rate \(\alpha\) controls the step size. If the rate is too large, the cost may oscillate or diverge. A rate that is too small may require many iterations. Recording the cost at regular intervals helps check optimization progress; a decreasing cost alone does not verify the gradients.

Vectorized prediction and training

Place \(m\) examples in columns:

\[X=\begin{bmatrix}|&|&&|\\x^{(1)}&x^{(2)}&\cdots&x^{(m)}\\|&|&&|\end{bmatrix}\in\mathbb{R}^{n_x\times m}\]

Then all predictions are calculated at once:

\[Z=w^TX+b,\qquad A=\sigma(Z)\]

NumPy broadcasts the scalar bias across all examples. Matrix operations handle the examples together, while the training loop repeats parameter updates.

For more detail on this programming pattern, see NumPy Vectorization and Broadcasting for Neural Networks.

Complete NumPy implementation

Run the following blocks in order: later blocks use the functions defined earlier. They train unregularized logistic regression on eight two-feature examples, then report probabilities and training-set metrics. Each update uses all eight examples.

1. Import NumPy and define the sigmoid function

import numpy as np


def sigmoid(z):
    """Numerically stable sigmoid for a NumPy array."""
    z = np.asarray(z, dtype=np.float64)
    result = np.empty_like(z)

    positive = z >= 0
    result[positive] = 1.0 / (1.0 + np.exp(-z[positive]))

    exp_z = np.exp(z[~positive])
    result[~positive] = exp_z / (1.0 + exp_z)
    return result

The two branches use equivalent formulas whose exponential arguments are nonpositive, avoiding exponential overflow for finite inputs. Floating-point sigmoid outputs can still round to exactly 0 or 1, so the loss will be calculated from the scores themselves.

2. Initialize parameters and calculate gradients

Zero initialization works for this single logistic unit: at the start each example contributes \(x(0.5-y)\) to the weight gradient, which need not vanish. There are no interchangeable hidden units that must be initialized differently.

To compute BCE from a score, substitute sigmoid into the loss and simplify: \(\mathcal L(z,y)=\log(1+e^z)-yz\). For positive \(z\), factor \(e^z\) out of the logarithm. Combining the two signs gives

\[\mathcal L(z,y)=\max(z,0)-yz+\log(1+e^{-|z|}).\]

The exponential now has a nonpositive argument, and np.log1p(u) evaluates \(\log(1+u)\) accurately even when \(u\) is small. This is the expression used below. At \(z=-800,y=1\), it returns 800 without taking the logarithm of a rounded probability.

def initialize_parameters(n_features):
    w = np.zeros((n_features, 1), dtype=np.float64)
    b = 0.0
    return w, b


def forward_and_backward(X, Y, w, b):
    """Return cost and gradients for one full-batch step."""
    assert X.ndim == 2 and X.shape[1] > 0
    m = X.shape[1]
    assert Y.shape == (1, m), "Y must have shape (1, m)"
    assert w.shape == (X.shape[0], 1) and np.ndim(b) == 0
    assert np.all((Y == 0) | (Y == 1)), "labels must be 0 or 1"
    assert all(np.isfinite(v).all() for v in (X, Y, w, np.asarray(b)))

    Z = w.T @ X + b
    A = sigmoid(Z)

    # Evaluate BCE directly from logits, without clipping probabilities.
    cost = np.mean(
        np.maximum(Z, 0) - Z * Y + np.log1p(np.exp(-np.abs(Z)))
    )

    dZ = A - Y
    dw = (X @ dZ.T) / m
    db = float(np.sum(dZ) / m)

    assert dw.shape == w.shape
    assert np.isfinite(cost)
    return cost, dw, db

One call returns the current cost and its parameter gradients. The input assertions make the shape convention explicit before broadcasting can combine incompatible inputs. These are teaching checks; a reusable API would normally raise descriptive exceptions for invalid inputs.

3. Train the model and define evaluation helpers

def train(X, Y, learning_rate=0.1, iterations=2000):
    w, b = initialize_parameters(X.shape[0])
    history = []

    for step in range(iterations):
        cost, dw, db = forward_and_backward(X, Y, w, b)

        w -= learning_rate * dw
        b -= learning_rate * db

        if step % 200 == 0 or step == iterations - 1:
            history.append((step, cost))

    return w, b, history


def predict(X, w, b, threshold=0.5):
    probabilities = sigmoid(w.T @ X + b)
    labels = (probabilities >= threshold).astype(int)
    return probabilities, labels


def binary_confusion_matrix(Y, predictions):
    y = Y.ravel()
    p = predictions.ravel()

    tn = int(np.sum((y == 0) & (p == 0)))
    fp = int(np.sum((y == 0) & (p == 1)))
    fn = int(np.sum((y == 1) & (p == 0)))
    tp = int(np.sum((y == 1) & (p == 1)))
    return np.array([[tn, fp], [fn, tp]])

During training, the same vectorized forward-and-backward calculation is repeated. The history list records occasional cost values without storing every iteration.

4. Create the dataset, train, and print the results

X = np.array([
    [-2.0, -1.5, -1.0, -0.5, 0.5, 1.0, 1.5, 2.0],
    [-1.0, -0.5, -1.5,  0.0, 0.5, 1.5, 1.0, 2.0],
])

Y = np.array([[0, 0, 0, 0, 1, 1, 1, 1]])

w, b, history = train(X, Y)
probabilities, predictions = predict(X, w, b)

accuracy = np.mean(predictions == Y)
matrix = binary_confusion_matrix(Y, predictions)

print("Initial recorded cost:", round(history[0][1], 6))
print("Last recorded cost:", round(history[-1][1], 6))
print("Weights:", np.round(w.ravel(), 4))
print("Bias:", round(b, 4))
print("Probabilities:", np.array2string(
    np.round(probabilities.ravel(), 4), suppress_small=True
))
print("Predictions:", predictions.ravel())
print("Accuracy:", accuracy)
print("Confusion matrix [[TN, FP], [FN, TP]]:")
print(matrix)

Verified output and interpretation

Running the code produces values close to the following:

Initial recorded cost: 0.693147
Last recorded cost: 0.013882
Weights: [4.439  2.6201]
Bias: -0.5808
Probabilities: [0.     0.0002 0.0001 0.0573 0.9502 0.9996 0.9998 1.    ]
Predictions: [0 0 0 0 1 1 1 1]
Accuracy: 1.0
Confusion matrix [[TN, FP], [FN, TP]]:
[[4 0]
 [0 4]]

The initial cost is approximately 0.693 because every probability is 0.5, giving a loss of −log(0.5) for either label. The displayed 0.0000 and 1.0000 probabilities are rounded values in this run. Training separates this deliberately simple dataset: the confusion matrix contains four true negatives and four true positives. The last recorded cost is measured before the final parameter update. np.array2string(..., suppress_small=True) keeps the rounded probabilities in decimal notation for this display.

This perfect training result does not establish generalization. The dataset is tiny and used only to check the implementation. It is also linearly separable, with no regularization: scaling a separating score can drive BCE toward zero without a finite minimizing parameter vector. Training stops at the specified iteration count, so weights and loss can keep changing after accuracy reaches 1. A real experiment uses validation data for choices such as the threshold, followed by evaluation on untouched test data with metrics chosen for the application.

Regularization and model limits

A common extension penalizes large weights with \(J_\lambda=J+\frac{\lambda}{2}\sum_j w_j^2\), where \(\lambda>0\) controls the strength and the intercept is left unpenalized. With this convention, add \(\lambda w\) to dw; db is unchanged. Some implementations scale the penalty by the sample size, so compare their objective definitions before copying a numeric strength.

The penalty trades training fit against coefficient size. For this dataset, which contains both classes, positive L2 regularization prevents the unbounded separating solution and gives a finite optimum. It can also reduce sensitivity to limited data, but stronger regularization need not improve held-out performance. Choose its strength on validation data. The code above intentionally leaves this penalty out; Ridge, Lasso, and Elastic Net develops these choices further.

Regularization does not change which features enter the score. If the required boundary cannot be represented by those features, shrinking their weights cannot supply the missing interaction or curvature. Inspect representation and validation errors alongside optimization progress.

Accuracy and the confusion matrix

Accuracy is the fraction of correct labels, but it can be misleading on imbalanced data. A classifier that always predicts the majority class can have high accuracy while failing to detect the class that matters.

Precision is \(TP/(TP+FP)\), the fraction of predicted positives that are correct; recall is \(TP/(TP+FN)\), the fraction of actual positives detected. Both require a nonzero denominator. For example, \(TP=8,FP=2,FN=4\) gives precision 0.8 and recall about 0.667. Changing the threshold changes these counts. Exercise 1 compares them and explicitly prices the errors. Ranking measures such as ROC-AUC use the scores across thresholds; their interpretation and probability calibration are covered in Classification Metrics and Thresholds.

Feature scaling

Features on very different scales can make gradient descent inefficient. For example, if one feature ranges from zero to one and another ranges into the millions, the cost surface may be highly elongated. A single learning rate may then need to be small to keep updates stable in the steep directions, slowing progress in others.

A common standardization uses the training-set mean and standard deviation:

\[x’_j=\frac{x_j-\mu_j}{\sigma_j}\]

Here \(\mu_j\) and \(\sigma_j\) are calculated from the training rows only; fit them separately within each training fold during cross-validation. Apply those same values to validation, test, and future inputs. Recomputing them on test data changes the fitted preprocessing rule and uses information from the evaluation sample. For a constant training feature, \(\sigma_j=0\); drop it or use a defined convention such as a scale of one instead of dividing by zero.

Numerical stability

A direct sigmoid implementation can overflow when it calculates \(e^{-z}\) for a very large negative \(z\). The implementation above uses separate formulas for positive and negative values. The loss is computed directly from logits using the stable expression derived above and tested at extreme scores in Exercise 3; it never takes the logarithm of a rounded sigmoid output.

The library function scipy.special.expit computes sigmoid values; it does not compute binary cross-entropy. Taking logarithms of its output can still encounter rounded probabilities of 0 or 1. Framework operations such as BCEWithLogitsLoss instead evaluate the loss directly from logits with a numerically stable formulation.

Debugging checklist

• Confirm that \(X\) has shape (n_features, m) and \(Y\) has shape (1, m).
• Confirm that \(w\), \(dw\), and their updates all have shape (n_features, 1).
• Verify that labels contain only zero and one.
• Check for NaN or infinite values in the cost, gradients, and parameters.
• Record the cost periodically; it should generally decrease for an appropriate learning rate.
• Compare one vectorized prediction with a manual scalar calculation.
• Do not evaluate only on the same examples used for training.
• If performance is poor, inspect feature scaling and class balance before assuming that the code is wrong.

Logistic regression as a one-unit neural network

Logistic regression can be drawn as a neural network with one output unit and no hidden layer:

\[x\longrightarrow z=w^Tx+b\longrightarrow a=\sigma(z)=\hat{y}\]

Forward propagation computes \(a\), the loss measures the error, and backpropagation produces \(dw\) and \(db\). Deeper networks repeat these ideas across more units and layers. That is why logistic regression is such a useful bridge between classical linear models and neural-network training.

The multiclass output layer and its fused gradient are covered in Softmax and Multiclass Classification Explained.

Exercises

1. Threshold is a choice. Given predicted probabilities and labels, compute accuracy, precision, and recall at thresholds \(0.3\), \(0.5\), and \(0.7\). Among these three candidates, which threshold maximizes accuracy, and which would you pick if a false negative costs ten times a false positive?

Use the synthetic scores below as a validation sample. They are not guaranteed to be calibrated probabilities. If no examples are predicted positive, report precision as zero for this exercise.

import numpy as np
rng = np.random.default_rng(0)
y = (rng.uniform(size=2000) < 0.2).astype(int)
p = np.clip(rng.normal(0.35 + 0.3*y, 0.2), 0, 1)

You should get: three different operating points, and a chosen threshold that is not the accuracy-maximizing one.

import numpy as np
rng = np.random.default_rng(0)
y = (rng.uniform(size=2000) < 0.2).astype(int)          # positive-label probability 0.2
p = np.clip(rng.normal(0.35 + 0.3*y, 0.2), 0, 1)
for t in (0.3, 0.5, 0.7):
    pred = (p >= t).astype(int)
    tp = ((pred==1)&(y==1)).sum(); fp = ((pred==1)&(y==0)).sum()
    fn = ((pred==0)&(y==1)).sum()
    print(t, "acc", round((pred==y).mean(),3),
          "prec", round(tp/max(tp+fp,1),3), "rec", round(tp/max(tp+fn,1),3),
          "FP", fp, "FN", fn, "cost", fp + 10*fn)
# 0.3 acc 0.53 prec 0.302 rec 0.95 FP 919 FN 21 cost 1129
# 0.5 acc 0.772 prec 0.473 rec 0.759 FP 354 FN 101 cost 1364
# 0.7 acc 0.837 prec 0.696 rec 0.394 FP 72 FN 254 cost 2612

2. Why not squared error, numerically. Compare the gradient magnitude of BCE and of MSE-with-sigmoid at \(z=-6\) with true label \(y=1\) — the confidently-wrong case. Report both and state which gives a larger correction to the logit under the same learning rate for a gradient step on this example.

You should get: one gradient magnitude near 1 and one near zero, differing by more than two orders of magnitude.

import numpy as np
z, y = -6.0, 1.0
a = 1/(1+np.exp(-z))                       # 0.00247
bce_grad = a - y                           # dL/dz for BCE
mse_grad = 2*(a - y) * a*(1-a)             # dL/dz for MSE + sigmoid
print(round(a,5), round(bce_grad,5), round(mse_grad,8))
# 0.00247 -0.99753 -0.00492082

3. Numerical stability. Compute BCE naively as \(-y\log\sigma(z)-(1-y)\log(1-\sigma(z))\) at \(z=-800\) and \(z=800\). Then compute it with the stabilized form \(\max(z,0)-zy+\log(1+e^{-|z|})\). Compare.

You should get: nan or inf from one form and a finite number from the other.

import numpy as np
def naive(z, y):
    a = 1/(1+np.exp(-z))
    return -(y*np.log(a) + (1-y)*np.log(1-a))
def stable(z, y):
    return np.maximum(z,0) - z*y + np.log1p(np.exp(-np.abs(z)))
for z in (-800.0, 800.0):
    print(z, naive(z,1.0), stable(z,1.0))
# -800.0 inf 800.0
# 800.0 nan 0.0