Logistic Regression Gradient Derived From Scratch (Full Math + NumPy)

There is one line of algebra hiding inside almost every classifier you will ever train. It says the gradient equals X transpose, times the predictions minus the labels, divided by n. In this lesson we will derive that line from nothing, so you never have to trust it again.

• Xᵀ(p − y) / n

As we saw in lesson 40, linear regression has a clean closed form. Logistic regression does not, so we need the gradient. Today we derive it, check it on three data points, extend it to softmax, then code gradient descent and Newton's method, watch the weights explode on separable data, and compare with scikit-learn.

• Derive the gradient
• Check by hand
• Softmax, Newton, NumPy

Recall from lesson 10 that a binary label is a Bernoulli variable. If the model says the probability of class one is p, then the probability of the observed label y is p to the power y, times one minus p to the power one minus y. When y is one, only p survives. When y is zero, only one minus p survives.

• P(y | x) = pʸ (1 − p)¹⁻ʸ
• p = σ(wᵀx)

Assuming the examples are independent, the likelihood of the whole dataset is the product of these terms over all n samples. A product of many numbers below one shrinks toward zero, and products are awkward to differentiate. So we do what we always do in maximum likelihood, and take the logarithm.

• Independent samples multiply
• L = ∏ pᵢʸⁱ (1 − pᵢ)¹⁻ʸⁱ

The log turns the product into a sum, and brings each exponent down in front. We get the sum over i of y times log p, plus one minus y times log of one minus p. Maximizing this is the same as maximizing the likelihood, because the log never changes where the peak is.

• log L = Σ y log p + (1−y) log(1−p)
• Maximize this

Machine learning prefers to minimize, so we flip the sign and average over n. This gives the loss J, equal to minus one over n, times the sum of y log p plus one minus y log of one minus p. That is exactly binary cross-entropy, the same quantity from lesson 10. Training by maximum likelihood and minimizing cross-entropy are one and the same thing.

• J = −(1/n) Σ [y log p + (1−y) log(1−p)]
• Binary cross-entropy

Look at one sample. If the label is one, the loss is minus log p. At p equals one the loss is zero, at p equal to one half it is about 0.69, and as p approaches zero it shoots to infinity. Confident wrong answers are punished brutally, and that is exactly the behavior we want.

• y = 1: loss = −log p
• y = 0: loss = −log(1 − p)

The model first computes a score z, equal to w transpose x plus b. Then the sigmoid squashes it into a probability: sigma of z equals one over one plus e to the minus z. Large positive scores give probabilities near one, large negative scores near zero, and zero gives exactly one half.

• σ(z) = 1 / (1 + e⁻ᶻ)
• z = wᵀx + b

Now we derive its slope in full. Write sigma as one plus e to the minus z, all raised to the power of negative one. The outer power rule gives negative one times that bracket to the power of negative two. The inner derivative of one plus e to the minus z is negative e to the minus z. The two negatives cancel.

• σ = (1 + e⁻ᶻ)⁻¹
• Chain rule on the outside

So the derivative is e to the minus z over one plus e to the minus z, squared. Split that into two fractions: one over one plus e to the minus z, which is sigma, times e to the minus z over one plus e to the minus z. That second fraction equals one minus sigma, because one minus sigma is the numerator e to the minus z over the same denominator. So d sigma d z equals sigma times one minus sigma.

• dσ/dz = e⁻ᶻ / (1 + e⁻ᶻ)²
• = σ · (1 − σ)

Sanity check. At z equals zero, sigma is 0.5, so the slope is 0.5 times 0.5, which is 0.25, the steepest point of the curve. At z equals two, sigma is about 0.881, so the slope is 0.881 times 0.119, about 0.105. The curve flattens as it saturates. Remember that flattening, because it matters later.

• z = 0: σ = 0.5, slope 0.25
• z = 2: σ = 0.881, slope 0.105

As we saw in lesson 36, a chain rule problem is easiest as a graph. The weights w feed into the score z, z feeds into the probability p, and p feeds into the loss. To get the derivative of the loss with respect to w, we multiply the local derivatives along that path.

• w → z → p → loss
• Multiply local derivatives

Step one. For a single sample the loss is minus y log p, minus one minus y times log of one minus p. Differentiate with respect to p. The first term gives minus y over p. The second gives plus one minus y over one minus p, because the inner derivative of one minus p is negative one.

• ℓ = −y log p − (1−y) log(1−p)
• ∂ℓ/∂p = −y/p + (1−y)/(1−p)

Put both terms over the common denominator p times one minus p. The numerator becomes minus y times one minus p, plus one minus y times p. Expand it: minus y plus y p, plus p minus y p. The y p terms cancel, leaving p minus y. So the derivative of the loss with respect to p is p minus y, over p times one minus p.

• ∂ℓ/∂p = (p − y) / (p(1−p))

Step two is the sigmoid derivative we just proved: d p d z equals p times one minus p. Multiply it by step one. The p times one minus p in the denominator cancels exactly with the one in the numerator. We are left with something wonderfully simple: the derivative of the loss with respect to the score z is p minus y.

• ∂ℓ/∂z = (p−y)/(p(1−p)) · p(1−p)
• = p − y

Step three. The score z is w transpose x plus b, so its derivative with respect to weight j is just feature j, and with respect to the bias it is one. Multiplying the chain gives the derivative of the loss with respect to w j equal to p minus y, times x j. The error, times the input.

• z = wᵀx + b
• ∂z/∂wⱼ = xⱼ

Now average over all n samples and stack them. The data matrix X has one row per sample. The vector p minus y holds every error. Multiplying X transpose by that vector sums error times input for each feature at once. So the gradient is X transpose times p minus y, divided by n, and the bias gradient is the mean of p minus y.

• ∇w J = Xᵀ(p − y) / n
• ∂J/∂b = mean(p − y)

Notice that this matches the shape of the gradient for linear regression from lesson 40: X transpose times prediction minus target. The only difference is that the prediction passes through a sigmoid. This is not a coincidence. It happens whenever you pair an exponential-family likelihood with its natural link, and cross-entropy was built for the sigmoid.

• Same form: Xᵀ(prediction − target)
• Different prediction function

Time for numbers. Take three points with a single feature: x equals minus one with label zero, x equals one with label one, and x equals two with label one. Start with w and b both zero. Every score is zero, so every probability is exactly 0.5, and the loss is log 2, about 0.693.

• (x, y): (−1, 0), (1, 1), (2, 1)
• Start w = 0, b = 0

The errors p minus y are 0.5, minus 0.5, and minus 0.5. For w, multiply each error by its x: minus 0.5, minus 0.5, and minus 1.0. These add to minus 2, divided by three gives about minus 0.667. For the bias, the errors add to minus 0.5, divided by three gives about minus 0.167.

• p − y = 0.5, −0.5, −0.5
• ∇w = −2/3, ∇b = −1/6

Take one step with a learning rate of one. Subtracting the negative gradient gives w equals 0.667 and b equals 0.167. The new scores are minus 0.5, 0.833, and 1.5. The probabilities are 0.378, 0.697, and 0.818. The new loss averages 0.474, 0.361, and 0.201, which is about 0.345. One step nearly halved the loss.

• w = 0.667, b = 0.167 (lr = 1)
• Loss: 0.693 → 0.345

Pause on something subtle. These three points can be split perfectly by a threshold between minus one and one. The two classes do not overlap at all. That innocent fact will cause real trouble later in this lesson, so keep it in the back of your mind.

• This tiny set is separable
• Remember this for later

Now generalize to K classes. Each class gets its own score z k, and softmax turns scores into probabilities: p k equals e to the z k, divided by the sum of e to the z j. The label is one-hot, and the loss is minus the sum of y k times log p k. With K equals two, this collapses back to the sigmoid case.

• pₖ = e^{zₖ} / Σⱼ e^{zⱼ}
• Loss = −Σ yₖ log pₖ

We need the derivative of p k with respect to z j. Using the quotient rule, there are two cases. If j equals k, we get p k times one minus p k, exactly like the sigmoid. If j differs from k, we get minus p k times p j. Both cases are written compactly as p k times, delta k j minus p j.

• ∂pₖ/∂zⱼ = pₖ(δₖⱼ − pⱼ)

The loss derivative with respect to p k is minus y k over p k. Multiply by the Jacobian and sum over k. The p k cancels, leaving minus the sum of y k times delta k j minus p j. The first part picks out y j. The second part is p j times the sum of the y k, which is one since the label is one-hot. The result is p j minus y j. Same shape as before.

• ∂ℓ/∂zⱼ = −Σₖ (yₖ/pₖ) pₖ(δₖⱼ − pⱼ)
• = pⱼ − yⱼ

Stacking everything, P is the n by K matrix of predicted probabilities and Y is the one-hot label matrix. The gradient for the weight matrix W is X transpose times P minus Y, divided by n. It is the same line as the binary case, just with a matrix of errors instead of a vector. One formula covers every linear classifier trained with cross-entropy.

• ∇W = Xᵀ(P − Y) / n
• P, Y: n × K matrices

Time to code. Computing log of sigma directly breaks for large scores, because e to the z overflows, or p rounds to exactly zero and log gives minus infinity. Rewrite the loss algebraically as softplus of z minus y times z, where softplus is log of one plus e to the z. Compute softplus as max of zero and z, plus log of one plus e to the minus absolute z. In NumPy that is logaddexp of zero and z.

• J = mean(softplus(z) − y·z)
• np.logaddexp(0, z)

The loop is short. Compute the scores z as X times w, with a column of ones in X for the bias. Compute p with a sigmoid. Update w by subtracting the learning rate times X transpose times p minus y, divided by n. Track the loss each step. Five lines of NumPy, with no automatic differentiation anywhere.

• p = σ(Xw)
• w ← w − lr · Xᵀ(p − y)/n

As we saw in lesson 11, plain gradient descent uses one step size for every direction. If the loss is steep in one direction and flat in another, a big step overshoots the steep one while a small step crawls along the flat one. We would like a method that knows the curvature. That is what the Hessian provides.

• Curvature differs by direction
• Learning rate is a compromise

Differentiate the gradient once more. Each p changes with z at rate p times one minus p, which we proved earlier. Collect those into a diagonal matrix S, with entries s i equal to p i times one minus p i. The Hessian is then X transpose S X, divided by n. Since every s i is positive, this matrix is positive semidefinite, so the loss is convex.

• g = Xᵀ(p − y)/n
• H = Xᵀ S X / n

Newton's method fits a quadratic bowl to the loss at the current point and jumps straight to its bottom. The update is w minus H inverse times g. In code, never form the inverse. Solve the linear system H delta equals g with np.linalg.solve. Near the optimum Newton converges extremely fast, often in under ten iterations.

• w ← w − H⁻¹ g
• Solve, do not invert

Rewrite the Newton step and something lovely appears. Each iteration solves a weighted least squares problem, like lesson 40, with weights s and a working response z tilde equal to X w plus S inverse times y minus p. That is why the method is called iteratively reweighted least squares. It is the classic algorithm behind statistical packages for generalized linear models.

• w = (XᵀSX)⁻¹ XᵀS ž
• ž = Xw + S⁻¹(y − p)

On a small problem, gradient descent needs hundreds of iterations to reach a loss tolerance, while Newton needs five or six. But each Newton step builds a d by d matrix and solves a system, which costs roughly d cubed. For a handful of features Newton wins easily. For thousands of features, gradient methods take over.

• GD: hundreds of cheap steps
• Newton: few steps, d³ cost each

Now back to that warning. If a line splits the two classes perfectly, we can always lower the loss by scaling the weights up. A bigger weight pushes every probability closer to its correct label. The loss shrinks toward zero but never reaches it, so there is no finite minimum to find.

• A line splits the classes perfectly
• Loss can approach zero

Run gradient descent on our three points and the weight keeps climbing forever, slowly, roughly like the log of the iteration count. Newton is even more dramatic. As probabilities saturate, the s values shrink toward zero, the Hessian becomes nearly singular, and the solver can warn about a singular matrix. Remember the flat sigmoid tails from earlier.

• ‖w‖ grows like log(iterations)
• Hessian entries → 0

The cure is to penalize large weights. Add lambda over two times the squared norm of w to the loss, as in lesson 17. The gradient gains one extra term, lambda times w, and the Hessian gains lambda times the identity, which also keeps it invertible. Now the loss has a finite minimum, and the weights settle at a sensible value.

• J + (λ/2)‖w‖²
• Gradient adds λw

Finally, verify. Fit scikit-learn's LogisticRegression with no penalty on a non-separable dataset, and compare with our gradient descent and Newton solutions. The coefficients should match to many decimal places. One trap: scikit-learn applies an L2 penalty by default, with C equal to one. Forget that, and your numbers will differ for a perfectly good reason.

• penalty=None, or C very large
• Coefficients agree

Let us recap. The negative log-likelihood of a Bernoulli model is binary cross-entropy. The sigmoid has derivative sigma times one minus sigma. The chain rule makes everything cancel to p minus y, giving X transpose p minus y over n. Softmax gives the same form, Newton uses the Hessian X transpose S X, and separable data sends the weights to infinity unless you regularize.

• NLL = cross-entropy
• dσ/dz = σ(1−σ)
• ∇ = Xᵀ(p−y)/n

Here is a teaser for the next lesson. The quantity p minus y, the negative gradient of the loss with respect to the score, will return as the target that each new tree in gradient boosting tries to fit. See you there.

• Lesson 42: Gradient Boosting and XGBoost
• Subscribe for the series