Logistic Regression — Fitting a Probability, and Deciding With It
A classification predicts a label, not a number. The smallest version of the problem — and by far the most common — is binary classification, where the answer can only be yes or no: malignant or benign, spam or not, fraudulent or legitimate. Logistic regression is the algorithm that solves it, and it solves it by refusing to answer the question directly. Instead of committing to a label, it outputs a probability, and only then does a threshold turn that probability into a decision.
Despite the name, logistic regression is a classification algorithm. It is called a regression for historical reasons only, and the name has stuck.
In one line. Logistic regression squashes the linear model through the sigmoid, $f_{\vec{w},b}(\vec{x}) = \sigma\left(\vec{w}\cdot\vec{x} + b\right)$, reads the output as $P(y = 1 \mid \vec{x})$, and fits $\vec{w}, b$ by minimising the log loss rather than the squared error — because a squared error wrapped around a sigmoid is not convex. The gradient that comes out is, remarkably, the same shape as linear regression's: $\frac{1}{m}\sum \left(f_{\vec{w},b}\left(\vec{x}^{(i)}\right) - y^{(i)}\right)x_j^{(i)}$.
What a classification predicts
A regression predicts a number drawn from infinitely many possible outputs — a price, a temperature. A classification picks from a small, fixed set of labels. When that set has exactly two members the problem is called binary classification, and $y$ can take only two values.
By convention those values are written $0$ and $1$. The class labelled $1$ is the positive class and the class labelled $0$ is the negative class — but positive here carries no suggestion of good. It marks the presence of the thing being detected: a malignant tumour is the positive class precisely because it is the thing you are looking for.
The obvious first idea is to reuse linear regression: fit a straight line and cut it at $0.5$. It fails, for two reasons.
- The output is unbounded. A line runs off to $\pm\infty$, so it happily predicts $-0.4$ or $3.7$ — neither of which is a label, and neither of which is a probability.
- The fit is dragged by points that are already correct. Add one training example far out along the positive side and the least-squares line tilts to reduce its (large, squared) error, which slides the $0.5$ crossing sideways and misclassifies examples that used to be right. A confidently-correct point should contribute nothing; under squared error it contributes the most.
What is needed is a model whose output is trapped between $0$ and $1$, which saturates instead of running away, and which can be read as a probability. That model is the sigmoid.
Defining the model
The sigmoid function
The sigmoid (or logistic) function is the S-shaped curve that maps the whole real line into the open interval between $0$ and $1$. It is the piece that turns a linear score into a probability, and it is the reason the algorithm has the name it does.
Its argument is written $z$ and it is deliberately centred on the origin: $z$ ranges over both signs, and $z = 0$ is the point of symmetry, where the curve passes through exactly $\tfrac{1}{2}$. Large positive $z$ drives the output toward $1$; large negative $z$ drives it toward $0$; and it never reaches either.
Two properties are worth naming now, because both are used later.
Both fall out of the algebra. For the first, $g(-z) = \frac{1}{1 + e^{z}}$; multiply top and bottom by $e^{-z}$ to get $\frac{e^{-z}}{e^{-z} + 1}$, which is exactly $1 - \frac{1}{1 + e^{-z}}$. It says the model is even-handed: the probability it assigns to the negative class at $-z$ is the probability it assigns to the positive class at $+z$.
For the second, write $g = \left(1 + e^{-z}\right)^{-1}$ and differentiate by the chain rule: $g' = -\left(1 + e^{-z}\right)^{-2} \cdot \left(-e^{-z}\right) = \frac{e^{-z}}{\left(1 + e^{-z}\right)^{2}}$. Split that as $\frac{1}{1 + e^{-z}} \cdot \frac{e^{-z}}{1 + e^{-z}}$ — the first factor is $g$, and the second is $1 - g$ by the identity just proved. The derivative of the sigmoid is therefore expressible entirely in terms of its own output, which is why the gradients in §6.2 collapse so cleanly.
The linear part, and the substitution
The sigmoid needs a number to squash, and that number is supplied by the familiar linear model. Its output is given the name $z$:
This is the entire connection between the two algorithms. Linear regression stops here and reports $z$ as the answer. Logistic regression treats $z$ as an intermediate — a raw, unbounded score, positive when the evidence points at the positive class — and passes it through the sigmoid to bound it.
So the two functions are not competing shapes fitted to the same points. They are stages of one model: the line computes the score, the sigmoid converts the score to a probability.
The logistic regression model
Substituting $z = \vec{w}\cdot\vec{x} + b$ into the sigmoid gives the algorithm.
The parameters are exactly the ones linear regression had — a weight vector $\vec{w}$ and a bias $b$ — and they play exactly the same roles. What changed is only what happens to their output on the way out of the model.
Notation
| Symbol | Meaning |
|---|---|
| $\vec{x}$, $y$ | input feature vector, output label with $y \in \{0, 1\}$ |
| $\vec{x}^{(i)}$, $y^{(i)}$ | the input and label of the $i$-th training example |
| $m$, $n$ | number of training examples; number of features |
| $x_j^{(i)}$ | the value of feature $j$ in the $i$-th training example |
| $\vec{w}$, $b$ | the parameters: the weight vector and the bias |
| $z$ | the linear score $\vec{w}\cdot\vec{x} + b$, the input to the sigmoid |
| $g(z)$ | the sigmoid, $1/(1 + e^{-z})$; also written $\sigma(z)$ |
| $f_{\vec{w},b}(\vec{x})$ | the model output — a probability in $(0, 1)$, not a label |
| $\hat{y}$ | the decision — the label the model commits to, $0$ or $1$ |
| $L$, $J$ | the loss on one example; the cost averaged over all $m$ |
| $\alpha$, $\lambda$ | the learning rate; the regularization parameter |
Keep $f$ and $\hat{y}$ apart. In linear regression they were the same object — the model's output was the prediction. Here they are two different things: $f_{\vec{w},b}$ is a probability, a real number strictly between $0$ and $1$, and $\hat{y}$ is the label you get after applying a threshold to it. Everything in §5 is about fitting $f$; all of §4 is about extracting $\hat{y}$ from it.
Interpreting the output
The output is a probability
Because the sigmoid never reaches its bounds, it is literally impossible for the model to output exactly $0$ or exactly $1$. That is not a defect to be engineered away — it is the point. The output is not a label; it is a degree of belief.
Read $f_{\vec{w},b}(\vec{x})$ as: there is an $f_{\vec{w},b}(\vec{x})$ chance that $y = 1$. A tumour that scores $0.7$ is one the model thinks has a $70$ percent chance of being malignant. And since the two classes are the only possibilities, the rest of the belief must sit on the other one — there is a $1 - f_{\vec{w},b}(\vec{x})$ chance that $y = 0$.
This is why only one number needs to be modelled even though there are two classes: the second is whatever is left over.
The formal notation
Written out with everything it is conditioned on, the model output is the probability that $y = 1$, given the input $\vec{x}$, and parameterised by $\vec{w}$ and $b$. The semicolon separates the two: what is observed, and what was learned.
This is the notation research papers use, and it is worth being fluent in: $\vec{x}$ is to the left of the semicolon because it is data the probability is conditioned on; $\vec{w}$ and $b$ are to the right because they are not random at all — they are the knobs the fitting procedure sets.
Reading it backwards: the log-odds
Given a probability $P = f_{\vec{w},b}(\vec{x})$, the odds are $\frac{P}{1 - P}$ — the chance it happens against the chance it does not. Inverting the sigmoid shows what the linear score $z$ was all along:
Where it comes from. Start from $P = \frac{1}{1 + e^{-z}}$ and solve for $z$. Taking reciprocals, $\frac{1}{P} = 1 + e^{-z}$, so $e^{-z} = \frac{1 - P}{P}$, and hence $e^{z} = \frac{P}{1 - P}$. Take logarithms and $z$ is the log-odds.
This is the sentence that makes logistic regression make sense. The model is linear — not in the probability, but in the log-odds of the probability. A one-unit increase in feature $x_j$ adds exactly $w_j$ to the log-odds, no matter where you are in feature space. The S-shape in Figure 1 is not a strange choice of curve; it is what a straight line looks like after you convert log-odds back into probabilities. The inverse map here is called the logit, and it is the reason the sigmoid — rather than some other S-shape — is the canonical choice.
The decision boundary
The decision, and the threshold
Sooner or later a probability has to become an action. The tumour is biopsied or it is not; the email lands in the inbox or in the spam folder. Turning $f_{\vec{w},b}(\vec{x})$ into a committed label $\hat{y}$ requires one more ingredient the model does not supply: a threshold.
The boundary is $z = 0$
Now ask when each case fires. With the default threshold of $0.5$, the question "when is $f_{\vec{w},b}(\vec{x}) \ge 0.5$?" unwinds one layer at a time:
The middle step is the one that matters, and it is free: the sigmoid crosses $0.5$ exactly at $z = 0$, so asking about the probability is the same as asking about the sign of the score. The whole non-linearity drops away at the boundary.
| Condition on the score | Probability | Decision |
|---|---|---|
| $\vec{w}\cdot\vec{x} + b \ge 0$ | $f_{\vec{w},b}(\vec{x}) \ge 0.5$ | $\hat{y} = 1$ |
| $\vec{w}\cdot\vec{x} + b < 0$ | $f_{\vec{w},b}(\vec{x}) < 0.5$ | $\hat{y} = 0$ |
The boundary in feature space
The set of points where the model is exactly undecided — where $z = 0$ — is called the decision boundary:
That is the equation of a straight line (in two features), a plane in three, a hyperplane in general. So although the model's output is curved, the frontier it draws through feature space is flat — the sigmoid bends the probabilities, not the boundary.
Moving the threshold
Nothing forces the threshold to be $0.5$. It is a choice, and choosing it is a statement about which mistake you would rather make.
Consider tumour detection, where the positive class is malignant. There are two ways to be wrong, and they are not remotely equal. A false positive sends a healthy patient for an unnecessary biopsy: unpleasant, expensive, survivable. A false negative sends a patient with cancer home: potentially fatal. Those costs are wildly asymmetric, so the decision rule should be too.
Lowering the threshold — flagging a tumour whenever $f_{\vec{w},b}(\vec{x}) \ge 0.2$, say — makes the model quicker to raise the alarm. It catches more of the real malignancies (fewer false negatives) at the price of more false alarms (more false positives). That is a trade worth making here, and a terrible one in a spam filter, where a false positive means a real email silently disappears.
The threshold is not learned. Gradient descent never sees it. Fitting produces the probability; the threshold is bolted on afterwards and can be retuned at any time — after deployment, per user, per hospital — without retraining a thing. Keeping the model's job ("how likely is this?") separate from the policy question ("what shall we do about it?") is the whole reason logistic regression outputs a probability rather than a label.
Non-linear boundaries
The boundary is a straight line only because $z$ is linear in the features. Feed the model polynomial features instead — $x_1^2$, $x_2^2$, $x_1 x_2$, and so on — and the same argument runs unchanged: the boundary is still the set where $z = 0$, but that set is now a curve.
With $z = w_1 x_1^2 + w_2 x_2^2 - 1$, for instance, the boundary $z = 0$ is the circle $x_1^2 + x_2^2 = 1$: inside it the model says $\hat{y} = 0$, outside it says $\hat{y} = 1$. Higher-order features bend the frontier into ellipses, and beyond that into shapes with no name at all.
Everything else in this note — the loss, the cost, the gradients, the update rule — is untouched by this. The features are just whatever numbers you chose to put in $\vec{x}$.
The cost function
Why not the squared error
Linear regression scored a candidate fit with the squared error cost:
The obvious move is to keep this and simply let $f_{\vec{w},b}$ be the new, sigmoid-wrapped model. It does not work. With a linear $f$, that cost is a perfect convex bowl — one minimum, and gradient descent cannot fail to find it. Wrap a sigmoid around $f$ and the bowl buckles: the surface grows bumps, ridges and multiple local minima, and gradient descent gets stuck in whichever one it happened to start above.
Why the sigmoid flattens the squared error. Take one example with $y = 1$ and let $a = g(z)$ be the model's output. The squared-error loss is $\frac{1}{2}\left(a - y\right)^{2}$, and its derivative with respect to the score is, by the chain rule and the identity in Equation 2, $\frac{\partial L}{\partial z} = \left(a - y\right) \cdot a\left(1 - a\right)$.
Look at what that does when the model is confidently wrong — $a \approx 0$ while $y = 1$. The error $\left(a - y\right)$ is as large as it can be, but the factor $a(1 - a)$ has collapsed to nearly zero, and the product with it. The gradient vanishes exactly where the model is worst. The surface goes flat far from the truth instead of steepening, and it is that flatness — the loss saturating instead of growing — that lets the extra bumps and the second minimum appear. The fix is to choose a loss whose derivative does not carry that $a(1 - a)$ factor, and the log loss is precisely the one that cancels it (see §6.2).
The loss for one example
Split the per-example term out of the cost and give it a name. The loss $L$ measures the model against a single training example; the cost $J$ is the average of the loss over the training set. Linear regression's loss was the squared error; logistic regression needs a different one.
The requirement is easy to state: the loss should be near zero when the prediction is near the target, and it should blow up as the prediction approaches the wrong answer. The logarithm does both.
Each branch is the negative logarithm of the probability the model assigned to the correct answer. When $y = 1$ the correct answer has probability $f$; when $y = 0$ it has probability $1 - f$. So the loss is one idea, not two: be surprised in proportion to how much probability you withheld from the truth.
The threat of infinity is the point. As the model drives its output toward the wrong label, the loss does not merely get large — it diverges. A model that says "definitely benign" about a malignant tumour is charged an unbounded price. The algorithm is therefore incapable of settling anywhere that is confidently wrong about a training example, and it is exactly this steepness that repairs the vanishing gradient of §5.1.
The simplified loss
Two cases are awkward to differentiate and worse to vectorize. They can be folded into a single expression:
It looks worse than what it replaces, and it is in fact the same thing. The trick is that $y^{(i)}$ is not a general number — it is $0$ or $1$, and nothing else. So in each term it acts as an on-off switch, and one of the two terms is always annihilated.
| If $y^{(i)}$ is… | First term | Second term | What is left |
|---|---|---|---|
| $1$ | $-1 \cdot \log(f)$ | $-(1 - 1)\log(1 - f) = 0$ | $-\log(f)$ |
| $0$ | $-0 \cdot \log(f) = 0$ | $-(1 - 0)\log(1 - f)$ | $-\log(1 - f)$ |
Each row of the table reproduces one branch of Equation 12 exactly. The single-line form buys two things: it can be differentiated once instead of twice, and it can be evaluated on a whole training set with array arithmetic and no branching.
The cost
The cost is the loss averaged over the training set — the same relationship it had in linear regression:
Substituting the simplified loss and pulling the shared minus sign out to the front gives the cost function in the form it is almost always written:
Mind the minus sign. Both logarithms are negative — their arguments are probabilities, strictly between $0$ and $1$ — so the bracket is negative, and the leading minus makes $J$ positive, as a cost must be. The sign can be carried inside each term instead, giving $\frac{1}{m}\sum \left[-y \log(f) - (1 - y)\log(1 - f)\right]$; it is the same function. What is not the same, and is a common slip, is putting a minus on only one of the two terms.
Where it comes from: maximum likelihood
This cost is not a guess that happened to work. It falls out of the maximum likelihood principle: choose the parameters under which the data you actually observed are as likely as possible.
The model says a single example is a coin flip landing on $1$ with probability $f$. Because $y^{(i)}$ is $0$ or $1$, the probability the model assigns to the label that was actually observed can be written in one expression, using the same on-off trick as before:
Assuming the examples are independent, the probability of the whole training set is the product of these over all $i$. Maximising a product of $m$ tiny numbers is numerically hopeless and algebraically unpleasant, so take the logarithm — which is monotonic, so it moves nothing — turning the product into a sum. Maximising that sum is the same as minimising its negative, and the negative log-likelihood, divided by $m$, is precisely Equation 15.
Two consequences follow, and they are the reason this cost is used rather than any other steep, well-behaved function:
- It is convex. The cost has a single global minimum and no local traps — the property Figure 3 showed the squared error losing. Gradient descent on it cannot get stuck.
- Its gradient is beautiful. The next section shows the sigmoid's derivative cancelling the logarithm's, leaving the plainest expression imaginable.
Gradient descent for logistic regression
The update rule
Fitting means finding the $\vec{w}$ and $b$ that minimise $J$. The procedure is unchanged from linear regression: start somewhere, compute the slope of the cost, and step downhill.
The derivatives
Here is where the choice of loss pays off. Work on a single example, write $z = \vec{w}\cdot \vec{x} + b$ and $a = g(z)$ for the model's output, and take the chain rule in two hops: from the loss to $a$, and from $a$ to $z$.
That cancellation is the whole design. The $a(1 - a)$ that the sigmoid contributes — the factor that went to zero when the model was confidently wrong, flattening the squared error and stranding gradient descent — is exactly the denominator that the logarithm's derivative produces. They annihilate, and what survives is $a - y$: the plain error, the amount by which the predicted probability missed the label. No saturation, no vanishing signal. A model that is confidently wrong now has $\partial L / \partial z \approx \pm 1$, the largest push the loss can give.
The log loss is, in this precise sense, the loss the sigmoid was asking for.
The rest is bookkeeping. Since $z = \vec{w}\cdot\vec{x} + b$, we have $\partial z / \partial w_j = x_j$ and $\partial z / \partial b = 1$; multiply through and average over the training set.
The final algorithm
for it in range(iterations):
z = X @ w + b # (m,) the linear scores
f = 1.0 / (1.0 + np.exp(-z)) # (m,) sigmoid -> probabilities
err = f - y # (m,) the error, a - y
dw = (X.T @ err) / m # (n,) one gradient per feature
db = err.mean() # scalar
w -= alpha * dw # simultaneous update
b -= alpha * db
The computation graph
The same algorithm read as a graph. The forward pass pushes one example through the two stages of the model to a loss; the backward pass walks the same edges in reverse, multiplying local slopes.
A run of the algorithm
Six tumours, one feature (size in centimetres). Sizes $1, 2, 3$ are benign ($y = 0$); sizes $4, 5, 6$ are malignant ($y = 1$). Start from $w = b = 0$ and run the updates with $\alpha = 0.5$.
| Iteration | $w$ | $b$ | $J(w,b)$ | boundary $x^{\ast}$ |
|---|---|---|---|---|
| $0$ | $0$ | $0$ | $0.6931$ | — |
| $1$ | $0.3750$ | $0$ | $0.6476$ | — |
| $2$ | $0.1876$ | $-0.1347$ | $0.5980$ | $0.718$ |
| $5$ | $0.2988$ | $-0.3576$ | $0.5557$ | $1.197$ |
| $10$ | $0.3803$ | $-0.7204$ | $0.5009$ | $1.894$ |
| $50$ | $0.8827$ | $-2.6625$ | $0.2939$ | $3.016$ |
| $100$ | $1.2503$ | $-4.0337$ | $0.2122$ | $3.226$ |
| $1000$ | $3.2223$ | $-11.1043$ | $0.0642$ | $3.446$ |
| $5000$ | $5.6616$ | $-19.6737$ | $0.0193$ | $3.475$ |
Three things in that table are worth pausing on.
- The first step moves $w$ only. At $w = b = 0$ every score is $z = 0$, so every output is $f = 0.5$. The bias gradient is the mean error, $\frac{1}{6}\sum(0.5 - y^{(i)})$, and with three examples of each class the errors cancel exactly: $\partial J / \partial b = 0$. The weight gradient does not cancel, because it weights each error by $x^{(i)}$, and the malignant tumours are the large ones. The model learns the direction before it learns the offset.
- $w$ then briefly goes backwards — from $0.375$ down to $0.188$. With $b$ still at $0$, a positive $w$ pushes every score up, so the three benign tumours are now over-predicted, and the cheapest immediate repair is to pull $w$ back while $b$ starts its long slide negative. Gradient descent optimises the pair, not each parameter in turn.
- The cost never stops falling. By iteration $5000$ the weights are still growing and $J$ is still shrinking, with no minimum in sight.
On separable data, the minimum is at infinity. This dataset is perfectly separable — a boundary at $x^{\ast} = 3.5$ classifies every example correctly. Once that is achieved, scaling $\vec{w}$ and $b$ up by a common factor leaves the boundary where it is but pushes every score further from zero, driving each output closer to its label and the log loss closer to zero. So $J$ can always be reduced a little more by making the weights a little larger, and the infimum $J = 0$ is only reached in the limit $\lVert\vec{w}\rVert \to \infty$. The fitted probabilities harden toward $0$ and $1$ — the model becomes ever more confident about training points it already classifies correctly.
In practice this is stopped by capping the iterations, or — properly — by regularization (§7), which puts a price on large weights and so restores a finite minimum. It is the cleanest possible motivation for the next section.
Same shape, different model
Set the two algorithms side by side and the update rules are, character for character, the same.
| Linear regression | Logistic regression | |
|---|---|---|
| Model $f_{\vec{w},b}(\vec{x})$ | $\vec{w}\cdot\vec{x} + b$ | $g\left(\vec{w}\cdot\vec{x} + b\right)$ |
| Output | a number, unbounded | a probability in $(0,1)$ |
| Predicts | the value $\hat{y}$ | $P(y = 1 \mid \vec{x})$, then threshold |
| Loss $L$ | $\frac{1}{2}\left(f - y\right)^2$ | $-y\log(f) - (1-y)\log(1-f)$ |
| Cost $J$ | $\frac{1}{2m}\sum \left(f - y\right)^2$ | $-\frac{1}{m}\sum \left[y\log f + (1-y)\log(1-f)\right]$ |
| $\partial J / \partial w_j$ | $\frac{1}{m}\sum \left(f - y\right) x_j$ | $\frac{1}{m}\sum \left(f - y\right) x_j$ |
| $\partial J / \partial b$ | $\frac{1}{m}\sum \left(f - y\right)$ | $\frac{1}{m}\sum \left(f - y\right)$ |
They look identical and they are not. The gradient expressions agree symbol for symbol, but $f_{\vec{w},b}$ means something different in each column — a line on the left, a sigmoid on the right. Substitute the definitions and the two are entirely different functions of $\vec{w}$ and $b$. What the coincidence really says is that the error times the input is the universal gradient signal for this family of models: it is what backpropagation carries through a deep network, layer after layer, and meeting it twice here is the first sign of that pattern.
Regularized logistic regression
A model with many features — especially polynomial ones, which bend the boundary of §4.5 into ever more elaborate shapes — can contort itself to separate the training set perfectly and generalise terribly. And as §6.5 showed, even a well-behaved separable dataset drives the weights to infinity if nothing stops them. The cure is the one linear regression used, applied without modification.
The regularized cost
Add a penalty proportional to the squared size of the weights. The only change to Equation 15 is the extra term on the end.
The parameter $\lambda$ sets the exchange rate between the two goals. At $\lambda = 0$ the penalty vanishes and the weights are free to run away; make $\lambda$ enormous and the cheapest thing to do is set every weight to zero, leaving a model that ignores its inputs and predicts the same probability for everyone. Useful values sit between.
The modified derivative
Differentiating the new term is immediate: $\frac{\partial}{\partial w_j} \frac{\lambda}{2m} \sum_k w_k^2 = \frac{\lambda}{m} w_j$, since every $w_k$ with $k \ne j$ is a constant as far as $w_j$ is concerned, and the $2$ from the square cancels the $2$ in $2m$ — which is why the $2$ is there.
The complete algorithm
Regrouping the weight update collects the two appearances of $w_j$ and shows what the penalty is really doing:
Every iteration begins by multiplying the weight by a number slightly less than one — with $\alpha = 0.1$, $\lambda = 1$ and $m = 50$ that factor is $0.998$ — a steady leak toward zero that the data must actively push back against. A weight survives only if the training set keeps paying for it. This is identical to the regularized update for linear regression; once again, only the meaning of $f_{\vec{w},b}$ differs.
Where this goes next
Logistic regression is the hinge between regression and neural networks, and almost everything about it survives the transition intact.
- It is already a neuron. A linear map followed by a non-linearity, trained by gradient descent on a loss — that is the whole content of Figure 5, and it is exactly one unit of a neural network. Stack these units in layers and the model in Equation 4 becomes the last layer of a deep network, with the sigmoid still on the output and the log loss of Equation 15 still measuring it. Nothing that follows replaces logistic regression; it wraps it in more of itself.
- The error is the signal. The cancellation in Equation 18 that reduced $\partial L / \partial z$ to the bare $a - y$ is not a curiosity of this model. It is the seed that backpropagation pushes backwards through every layer of a deep network — the quantity called $dz$ in the reference — and the reason the sigmoid and the log loss are paired the way they are.
- Beyond two classes. With $N$ classes the sigmoid generalises to the softmax and the binary cross-entropy of Equation 13 to the general cross-entropy. The structure — score, squash, negative log-likelihood of the true class — does not change at all, and neither does the gradient: $a - y$ becomes the vector $\vec{p} - \vec{y}$. That generalisation is worked out in full in multiclass-softmax.md, which derives the softmax from the same log-odds argument used here, and shows that setting $N = 2$ collapses it back to this note.
The companion note neural-networks-reference.md picks the story up at exactly this point: it opens with logistic regression as the single neuron, and follows it to shallow networks, deep networks, and the component-level view of backpropagation. The note on linear regression is the story underneath this one.