Logistic Regression and Softmax: Derived from Scratch, Step by Step

Every time your email app decides whether a message is spam, it answers a yes or no question in a fraction of a second. Behind that decision, in many systems, is a surprisingly simple idea: a straight line, squashed through an S shaped curve. Today we derive that idea completely, so nothing in it feels like magic.

• Yes or no decisions
• Line plus squash
• Derived from scratch

Last episode, we built linear regression and fitted it with both the normal equation and gradient descent. There, the model predicted a number, like a house price. In this episode the model must predict a category, and that small change forces us to rethink the output, the loss, and the gradient.

• Last time: predict numbers
• Today: predict categories

Here is the plan. We start with classification versus regression, then meet odds and log-odds, and derive the sigmoid function from them. Next we study the decision boundary, derive binary cross-entropy from maximum likelihood, and prove the gradient step by step. Finally we code it, compare with scikit-learn, and extend everything to many classes with softmax.

• Odds to sigmoid
• Cross-entropy and gradient
• Code, then softmax

Regression predicts a continuous quantity, such as a temperature or a price. Classification predicts a discrete label, such as spam or not spam, or cat, dog, or bird. For two classes we encode the label y as either zero or one, and our goal becomes estimating the probability that y equals one, given the features x.

• Regression: continuous output
• Classification: discrete label
• Label encoded as 0 or 1

Suppose students study for some hours, and y is one if they pass. Fit a straight line by least squares and look at the result. For very long study times the line predicts a value above one, and for very short times it goes below zero. Those are not probabilities, so the output of a plain line cannot be trusted as one.

• Line goes above 1
• Line goes below 0
• Not a valid probability

There is a second problem. Add one student who studied for fifty hours and passed. That point is obviously a pass, yet it tilts the least squares line, and the place where the line crosses zero point five slides to the right. A single easy example should not change our decisions about the hard cases, so we need a different kind of model.

• One far point tilts the line
• Decision cutoff drifts

So the requirements are clear. We want an output always between zero and one, we want it to rise smoothly as evidence for class one grows, and we want it built from a linear score so that training stays simple. To find the right squashing function, we will not guess. We will derive it, starting from the idea of odds.

• Output between 0 and 1
• Smooth, monotone rise
• Linear score inside

Odds compare the chance something happens to the chance it does not. If p is the probability, the odds are p divided by one minus p. A probability of zero point eight gives odds of zero point eight over zero point two, which is four. We say four to one in favor.

• odds = p / (1 - p)
• p = 0.8 gives odds 4

Probability zero point five gives odds of one, a perfectly even bet. Probability zero point two gives odds of zero point two five. Probability zero point nine gives odds of nine. Notice that odds live between zero and infinity, which is better than the zero to one range, but they are lopsided: everything in favor of the event is crammed between zero and one.

• p = 0.5 gives odds 1
• p = 0.2 gives odds 0.25
• p = 0.9 gives odds 9

To fix the lopsidedness, take the natural logarithm. The log-odds, also called the logit, equal the natural log of p over one minus p. Now odds of four become about one point three eight six, odds of one become exactly zero, and odds of one quarter become negative one point three eight six.

• logit(p) = ln(p / (1 - p))
• Odds 4 gives 1.386
• Odds 1/4 gives -1.386

Log-odds are symmetric: swapping the two classes just flips the sign. And they stretch across the entire real line, from negative infinity to positive infinity, exactly like the output of a linear function. As probability goes to zero, the logit heads to minus infinity. As probability goes to one, it heads to plus infinity.

• Symmetric around zero
• Range: all real numbers
• Matches a linear score

Here is the one assumption that defines logistic regression. We assume the log-odds of class one are a linear function of the features. In symbols, the log of p over one minus p equals z, where z equals w transpose x plus b. We are not claiming the probability is linear, only the log-odds.

• Log-odds = w transpose x plus b
• Linear in log-odds, not in p

Now we invert the logit. Start with the natural log of p over one minus p equals z. Exponentiate both sides, so p over one minus p equals e to the z. Multiply out: p equals e to the z times one minus p, which is e to the z minus e to the z times p.

• Exponentiate both sides
• Expand the right-hand side

Move the term with p to the left: p plus e to the z times p equals e to the z. Factor out p, giving p times one plus e to the z equals e to the z. Divide, and p equals e to the z over one plus e to the z. Divide top and bottom by e to the z, and we get one over one plus e to the minus z. That is the sigmoid.

• Factor out p
• Divide by e to the z
• Result: 1 / (1 + e^-z)

Let us plug in numbers. At z equals zero, the sigmoid is exactly one half. At z equals two, it is one over one plus e to the minus two, about zero point eight eight one. At z equals minus two it is about zero point one one nine, and the two always add up to one. Large positive z gives values near one, large negative z near zero.

• sigma(0) = 0.5
• sigma(2) is about 0.881
• sigma(-2) is about 0.119

The sigmoid has a derivative we will need soon. Write it as one plus e to the minus z, all raised to the power of minus one. The chain rule gives e to the minus z over one plus e to the minus z, squared. Split that into sigma times e to the minus z over one plus e to the minus z, and the second factor equals one minus sigma. So the derivative is sigma times one minus sigma.

• sigma prime = sigma(1 - sigma)
• Maximum slope 0.25 at z = 0
• Flat when saturated

The complete model has two stages. First compute the linear score z equals w transpose x plus b. Then pass it through the sigmoid to get the probability of class one. Try one feature, hours studied, with weight one point five and bias minus three. At three hours, z is one point five and the probability is about zero point eight one eight.

• z = w x + b
• p = sigma(z)
• x = 3 gives p = 0.818

To turn a probability into a label, we pick a threshold, commonly zero point five. The sigmoid equals one half exactly when z equals zero. So the model predicts class one whenever w transpose x plus b is positive, and class zero when it is negative. The decision boundary is the set of points where z is zero.

• Threshold 0.5 means z = 0
• z > 0 predicts class 1
• Boundary: w x + b = 0

Back to our student example with weight one point five and bias minus three. Setting one point five x minus three equal to zero gives x equals two. Two hours is the boundary. At one hour the probability is about zero point one eight two, at two hours exactly one half, and at three hours about zero point eight one eight.

• 1.5x - 3 = 0 gives x = 2
• x = 1: p = 0.182
• x = 3: p = 0.818

With two features, the boundary w one x one plus w two x two plus b equals zero is a straight line. Take w equal to one, one, and b equal to minus three, so the boundary is x one plus x two equals three. The weight vector points perpendicular to this line, toward the class one side. The point two, two gives z equals one and probability about zero point seven three one.

• Boundary is a line
• w is perpendicular to it
• (2,2): z = 1, p = 0.731

The size of the weights controls how sharply the probability changes. Doubling w and b keeps the same boundary but makes the sigmoid steeper, so the model becomes more confident. Also note that this boundary is always flat: a line in two dimensions, a plane in three. If classes form a ring, the model needs extra features such as squared terms, which we meet in a later episode.

• Larger weights, steeper curve
• Same boundary, more confidence
• Boundary is always linear

We now need a loss function. Instead of guessing, we use the principle from our maximum likelihood episode: choose parameters that make the observed data as probable as possible. Each label is a coin flip, a Bernoulli variable, whose chance of heads is the model's probability p. So we simply ask how likely our labels are under the model.

• Labels are Bernoulli coin flips
• Maximize data probability
• Recall: maximum likelihood

The probability of one label is p if y equals one, and one minus p if y equals zero. We can write both cases in a single expression: p to the power y, times one minus p to the power one minus y. Check it: when y is one, the second factor has exponent zero and disappears, leaving p. When y is zero, the first factor disappears, leaving one minus p.

• P(y) = p^y (1-p)^(1-y)
• y = 1 leaves p
• y = 0 leaves 1 - p

With n independent examples, the likelihood is the product of these terms. Products of small numbers are awkward, so take the logarithm, which turns the product into a sum: the sum of y log p plus one minus y times log of one minus p. Maximizing this is the same as minimizing its negative. Divide by n to get the average, and we have binary cross-entropy.

• Likelihood is a product
• Log turns it into a sum
• Negate and average: BCE

Look at the loss for a single example whose true label is one. If the model says zero point eight, the loss is minus the natural log of zero point eight, about zero point two two three. At zero point one, it is about two point three zero three. At zero point zero one, it is about four point six zero five. Confident wrong answers are punished heavily, and the loss grows without bound as the probability approaches zero.

• p = 0.8: loss 0.223
• p = 0.1: loss 2.303
• p = 0.01: loss 4.605

A common mistake is to combine the sigmoid with mean squared error. The result is a non-convex loss surface with flat regions. When the model is confidently wrong, the sigmoid is saturated, its slope is nearly zero, and the gradient almost vanishes, so learning stalls. Cross-entropy cancels that slope, as we are about to see, and for logistic regression it is convex, so gradient descent has a single global minimum.

• MSE plus sigmoid: flat regions
• Cross-entropy is convex
• Single global minimum

Now for the central derivation. For one example, the loss L depends on p, p depends on z, and z depends on w. So we apply the chain rule through three links: the derivative of L with respect to p, times the derivative of p with respect to z, times the derivative of z with respect to w. We will compute each link separately.

• L depends on p
• p depends on z
• z depends on w

For one example, L equals minus y log p minus one minus y times log of one minus p. The derivative of log p is one over p, and the derivative of log of one minus p is minus one over one minus p. So dL by dp equals minus y over p, plus one minus y over one minus p. Combine over a common denominator, and it becomes p minus y, divided by p times one minus p.

• d ln p = 1 / p
• Common denominator
• dL/dp = (p - y) / (p(1-p))

The second link is the sigmoid derivative we proved earlier. The derivative of p with respect to z is p times one minus p. Notice something wonderful: this is exactly the denominator we just got in the first link, so the two will cancel.

• dp/dz = p(1 - p)
• Same factor as before
• They cancel

The third link is simple: z equals w transpose x plus b, so the derivative of z with respect to w is just x, and with respect to b is one. Multiplying the three links, the p times one minus p factors vanish, leaving dL by dw equals p minus y, times x. And dL by db equals p minus y.

• dz/dw = x, dz/db = 1
• dL/dw = (p - y) x
• dL/db = p - y

Interpret it: p minus y is simply the prediction error. If the label is one and the model says zero point two, the error is negative zero point eight, and gradient descent pushes the weights to raise the score along that example's features. This is the same form as linear regression, prediction minus target, times input. The sigmoid and cross-entropy fit together perfectly.

• p - y is the error
• Same form as linear regression
• Sigmoid slope cancelled

Take one example with x equal to two and y equal to one, and start with w and b both zero. Then z is zero, p is one half, and the loss is about zero point six nine three. The gradient for w is minus zero point five times two, which is minus one, and for b it is minus zero point five. With learning rate zero point one, w becomes zero point one and b becomes zero point zero five. The new z is zero point two five, p rises to about zero point five six two, and the loss falls to about zero point five seven six.

• Start: p = 0.5, loss 0.693
• Gradient: -1 and -0.5
• After step: loss 0.576

For a whole dataset, stack the inputs into a matrix X with n rows. Compute all scores at once, apply the sigmoid to get a vector p, and the average gradient is one over n times X transpose times the vector p minus y. The bias gradient is the mean of p minus y. Just a few lines of NumPy, no loops over examples.

• p = sigma(Xw + b)
• grad w = X^T (p - y) / n
• grad b = mean(p - y)

Here is the implementation in words. Initialize the weights to zeros and the bias to zero. In each iteration, compute z as X times w plus b, apply the sigmoid, compute the error p minus y, then form the gradients with X transpose times the error divided by n. Update the weights by subtracting the learning rate times the gradient, and do the same for the bias.

• Init w and b to zero
• Forward, error, gradient
• Update with learning rate

Run this on a two feature dataset with two clusters, and track the loss every iteration. It falls quickly at first, then flattens as it approaches the minimum. Meanwhile the decision line rotates and shifts until it cuts cleanly between the clusters. Plotting the loss curve is the easiest way to confirm the gradient is correct and the learning rate is sensible.

• Loss falls then flattens
• Boundary rotates into place
• Plot the loss curve always

First, numerical safety: if p becomes exactly zero or one, the logarithm blows up, so clip probabilities by a tiny epsilon when computing the loss. Second, scale your features, otherwise the loss surface is a stretched valley and gradient descent zigzags. Third, if the classes are perfectly separable, the weights grow without limit, because more confidence always lowers the loss. Regularization, coming in episode seventeen, solves that.

• Clip log inputs
• Scale the features
• Separable data: weights explode

Now fit scikit-learn's LogisticRegression on the same data. One detail matters: by default it applies L2 regularization with strength C equal to one. To compare fairly, turn the penalty off or set C very large. Then the learned coefficients and intercept match our from scratch values closely, differing only by tiny optimizer tolerances, since scikit-learn uses a different solver than plain gradient descent.

• Default has L2, C = 1
• Disable penalty to compare
• Coefficients nearly identical

Remember that the model's real output is a probability. In scikit-learn, predict returns labels using the zero point five threshold, while predict_proba returns the probabilities themselves. Keep them: they let you rank cases, choose a different threshold, and measure how well calibrated the model is. We will study exactly those topics in the next episode.

• predict gives labels
• predict_proba gives probabilities
• Threshold is a choice

What if there are three classes, or a hundred? A single sigmoid only handles two. The natural extension is to compute one score, or logit, for each class, using one weight vector per class. We stack them in a matrix W with K rows, so the score vector is z equals W x plus b. Now we need a function that turns K scores into K probabilities.

• One score per class
• z = W x + b
• Need K probabilities

Softmax does that. The probability of class k equals e to the z k, divided by the sum of e to the z j over all classes. Take scores two, one, and zero. The exponentials are about seven point three eight nine, two point seven one eight, and one, which sum to eleven point one oh seven. The probabilities are about zero point six six five, zero point two four five, and zero point zero nine zero.

• p_k = e^(z_k) / sum e^(z_j)
• Scores 2, 1, 0
• Probabilities 0.665, 0.245, 0.090

The exponential makes every number positive, the division makes them sum to one, and the result keeps the order of the scores: a bigger score always gets a bigger probability. Softmax is a soft version of taking the maximum, since a large gap between scores pushes the winner toward probability one. Adding the same constant to every score changes nothing, because the constant cancels in the ratio.

• Positive and sums to 1
• Soft version of argmax
• Shift-invariant

That shift invariance rescues us from overflow. With scores of one thousand, one thousand and one, and one thousand and two, computing the exponential directly overflows to infinity in floating point. Instead, subtract the largest score first, giving minus two, minus one, and zero. The probabilities come out the same, about zero point zero nine, zero point two four five, and zero point six six five, with no overflow.

• exp(1000) overflows
• Subtract the max score
• Same probabilities, stable

With only two classes, softmax collapses into the sigmoid. The probability of class one is e to the z one over e to the z one plus e to the z two. Divide top and bottom by e to the z one, and it becomes one over one plus e to the minus of z one minus z two. So it is a sigmoid of the score difference, and logistic regression is just the two class special case.

• Two classes: softmax = sigmoid
• Depends on z1 - z2
• Logistic regression is the special case

The loss generalizes naturally. Encode the true class as a one-hot vector, for example one, zero, zero. The loss is minus the sum over classes of y k times log p k, which keeps only the log probability of the true class. If the true class is zero and p is zero point six six five, the loss is about zero point four zero eight. This is exactly maximum likelihood for a categorical distribution.

• Target is one-hot
• L = -sum y_k ln p_k
• Only true class counts

Differentiate this loss with respect to the scores, and the same miracle appears: the gradient with respect to z k equals p k minus y k. For our example that is minus zero point three three five, plus zero point two four five, and plus zero point zero nine zero. The score of the true class is pushed up, and the others are pushed down. The weight gradient is that error vector times x transposed, the exact analogue of the binary case.

• dL/dz_k = p_k - y_k
• True class pushed up
• Others pushed down

There is a second way to handle many classes. In one-versus-rest, train K separate binary logistic regressions, where model k asks whether an example is class k or anything else. At prediction time, run all K models and pick the class with the highest score. This reuses everything we derived for two classes.

• Train K binary models
• Class k versus the rest
• Pick the highest score

What is the difference? One-versus-rest trains each model independently, so their probabilities do not sum to one and need rescaling. Each model also faces an imbalanced problem, one class against everything else. Softmax regression trains all classes jointly with a single loss and gives properly normalized probabilities. In practice softmax is usually preferred for logistic style models, while one-versus-rest is handy when you only have a binary algorithm to build on.

• OvR: independent, not normalized
• Softmax: joint, sums to 1
• Softmax usually preferred

Let us recap the whole chain. Classification predicts categories, and a plain line fails because it is not a probability. Assuming linear log-odds gives the sigmoid, with a linear decision boundary where z is zero. Maximum likelihood gives binary cross-entropy, whose gradient is simply p minus y, times x. Softmax extends this to many classes with the same elegant gradient, and one-versus-rest offers an alternative.

• Log-odds linear gives sigmoid
• MLE gives cross-entropy
• Gradient is (p - y) x
• Softmax for many classes

But a trained classifier raises new questions. Is ninety percent accuracy actually good? What threshold should we choose, and can we trust the probabilities? In the next episode, we cover confusion matrices, precision, recall, ROC curves, and calibration.

• Next: metrics and calibration
• Episode 16
• Subscribe for more