Logistic Regression,
built from scratch.

01 Starting Point: Connection with Linear Regression

Logistic regression is easy to mistake for a small variation on linear regression, and honestly, for the first half of the pipeline, that's true. Both start by combining the input features into a single linear score:

Everything up to computing z is identical: the same features X, the same kind of weights w, the same bias b, the same matrix multiplication. If you've implemented linear regression, you already know how to build this part of logistic regression.

The two models diverge right after z is computed. Linear regression treats z itself as the prediction. Logistic regression refuses to do that, because z is an unbounded real number — it can be -500 or 3000 — and that's not a usable answer to a yes/no question like "will this person buy a house?" A prediction for a binary outcome needs to live on a scale that means something, like a probability between 0 and 1. So logistic regression pushes z through one more function — the sigmoid — before calling anything a "prediction."

That one extra step changes everything downstream: the loss function, the interpretation of the output, and the shape of the gradient.

02 The Sigmoid Function

Sigmoid is the function that turns the raw linear score into something that can be read as a probability:

Why this particular function? Because it takes any real number and maps it into the open interval (0, 1), and it does so smoothly and monotonically — bigger z always means bigger (or equal) p, never a reversal. That monotonic squashing is exactly what lets p be interpreted as \(P(y=1 \mid x)\): the model's belief that the positive class is correct, given the input.

Three anchor points are worth memorizing, and the widget above lets you feel them directly:

• Large negative z → e^{-z} explodes → the fraction shrinks toward 0 → p → 0 (confident "no").
• z = 0 → e^{0} = 1 → p = 1/2 = 0.5 (total uncertainty).
• Large positive z → e^{-z} → 0 → p → 1 (confident "yes").

03 Where Log Loss Actually Comes From: Bernoulli → Likelihood

It's easy to treat -[y log(p) + (1-y) log(1-p)] as a formula that was simply handed down. It wasn't. It falls directly out of modeling each label as a coin flip whose bias is p — a Bernoulli distribution:

Check it against the two possible labels: if y = 1, the exponents collapse the expression to p. If y = 0, it collapses to 1-p. One formula, both cases — that's the trick of writing it with exponents on y and 1-y.

Multiplying these probabilities across a whole dataset gives the likelihood of the data under the current weights. Products of many small numbers between 0 and 1 shrink toward zero and are numerically miserable to work with, so we take the log to turn the product into a sum:

Training means finding the weights that make the observed labels as likely as possible — maximum likelihood estimation. Optimizers are built to minimize things, not maximize them, so we flip the sign and minimize the negative log likelihood instead. That flipped, per-example quantity is exactly log loss (binary cross-entropy):

This form also explains something that feels intuitive but is worth deriving: confidently wrong predictions are punished hard. Compare:

• y = 1, p = 0.9 → loss = -log(0.9) ≈ 0.105 — small, the model was right and confident.
• y = 1, p = 0.1 → loss = -log(0.1) ≈ 2.303 — over 20× larger, because the model was confidently wrong.

The logarithm is what creates that asymmetry — as p approaches 0 while the true label is 1, -log(p) shoots toward infinity. This is by design: it forces the model to be humble when it isn't sure.

04 What Comes From the Dataset vs. What Comes From the Model

A recurring source of confusion while debugging was blurring together things that come from the data and things the model computes. Keeping this boundary sharp makes every later formula easier to read.

p is the model's guess at the answer; y is the dataset's actual answer. The entire point of the loss function is to compare these two things — a computed guess against a given truth — and the entire point of training is to nudge w and b until the guesses land closer to the truth.

05 Matrix Shapes — and the Orientation I Kept Getting Backwards

With m training examples and n features, the shapes that make every later equation dimensionally consistent are:

Building X from separate feature arrays with np.column_stack is where the orientation clicked for me:

age = [20, 30, 40]
income = [20000, 30000, 40000]
data = np.column_stack((age, income))

So data[i][j] reads as: example i, feature j. The widget below makes the two axes physically obvious — click a row button to select one full training example, or a column button to select one feature across every example.

06 Loss vs. Gradient — Two Different Questions

These two quantities get used together so often that it's tempting to think of them as one thing. They aren't:

• Loss — \(J(w,b)\) — answers "how bad are my current predictions?" It's a single number summarizing performance right now.
• Gradient — \(\nabla J\) — answers "which direction should I change the parameters to make the loss smaller?" It's a direction, not a summary.

The cleanest analogy is an ordinary function like \(f(x) = x^2\). The value f(x) tells you your current height on the curve. The derivative f'(x) = 2x tells you which way is downhill from where you're standing. Loss is the height; gradient is the slope under your feet.

07 Deriving the Gradient, One Chain-Rule Link at a Time

Rather than accept dw = (p - y) * x on faith, it's worth walking the chain that produces it. The weight w_j influences the loss only by first influencing z, which influences p, which finally influences L:

Three functions are stacked here:

Because w_j only reaches L through this chain, the chain rule says its overall effect is the product of each link's local effect:

Each piece is manageable on its own:

1. \(\dfrac{\partial L}{\partial p}\) — differentiating -[y log p + (1-y) log(1-p)] with respect to p gives -y/p + (1-y)/(1-p).

2. \(\dfrac{\partial p}{\partial z}\) — the derivative of sigmoid has a famously tidy closed form:

3. \(\dfrac{\partial z}{\partial w_j}\) — since z = w_1x_1 + w_2x_2 + \dots + b, the only term containing w_j is w_j x_j, so this derivative is simply x_j.

Multiplying all three together, the p(1-p) terms cancel beautifully against the denominators from step 1:

$$\left(\frac{-y}{p} + \frac{1-y}{1-p}\right)\cdot p(1-p)\cdot x_j \;=\; (p - y)\,x_j$$

That cancellation is the whole reason logistic regression's gradient looks so clean — sigmoid and log loss were built to fit together like this.

Averaged over the whole dataset, in matrix form:

The term (p - y) — the model's probability minus the true label — is the model's error. It appears naturally, not by design choice: it's what's left after the calculus, and it has an intuitive reading too — if the model is perfectly calibrated, error is zero and the gradient vanishes; the more wrong p is, the larger the push to correct w.

08 Why the \(\frac{1}{m}\)?

The overall loss used for training is the average of the per-example losses, not their raw sum:

Differentiating an average just differentiates each term and keeps the same averaging constant out front — the 1/m doesn't depend on w, so it passes straight through:

Substituting the single-example gradient derived above gives the familiar batch formula:

Two counts show up in this codebase and they are not interchangeable:

• m — number of training examples (rows of X) — this is what the gradient divides by.
• len(param) — number of features / parameters — this has nothing to do with averaging the gradient.

Technically, nothing stops you from using the raw sum of losses instead of the mean. But then the scale of the gradient would depend on how many training examples you happened to have — a dataset with 10,000 rows would produce gradients roughly 10× larger than the same problem with 1,000 rows, for no reason related to the actual difficulty of the problem. Dividing by m keeps the gradient's scale — and therefore the learning rate that works well — independent of dataset size.

09 Gradient by Loops — the Literal Version

Before trusting a vectorized one-liner, it's worth writing the gradient exactly as the math reads: one accumulator per parameter, summed across every example.

error = p - y

dw = []

for j in range(len(param)):
    gradient = 0

    for i in range(len(data)):
        gradient += error[i] * data[i][j]

    gradient /= len(data)
    dw.append(gradient)

db = np.sum(error) / len(data)

The inner line is the entire derivation from Section 7, spelled out per index:

Summing that product over every example i, then dividing by len(data), produces exactly one number: the gradient for parameter j. Doing that for every j fills out the full dw vector.

Two indexing mistakes are easy to make here and both come from the row/column confusion in Section 5:

• data[j] when j is meant to be a feature index is wrong — data[j] accesses row j (a whole training example), not a feature column.
• data[:j] is wrong for a different reason — it's a slice from the start up to row j (still selecting rows), not a single column at all.

The only correct way to grab a full feature column is data[:, j] — keep every row (:), fix the column (j).

10 The Vectorized Gradient

Once the loop version is trusted, NumPy can do the same accumulation in one shot using matrix multiplication:

error = p - y

dw = (data.T @ error) / len(data)

db = np.sum(error) / len(data)

Matching the shapes makes it clear why this reproduces the loop exactly:

Transposing data turns "rows = examples" into "rows = features," so each row of data.T is one feature's values across all m examples. Multiplying that row by the length-m error vector and summing (which is exactly what matrix multiplication does) reproduces the inner loop's sum over i — but for every feature j simultaneously. The division by len(data) at the end is the same averaging step from Section 8.

error = p - y
dw = []
for j in range(len(param)):
    gradient = 0
    for i in range(len(data)):
        gradient += error[i] * data[i][j]
    gradient /= len(data)
    dw.append(gradient)
db = np.sum(error) / len(data)

error = p - y
dw = (data.T @ error) / len(data)
db = np.sum(error) / len(data)

11 Gradient Descent — Actually Using the Gradient

The gradient only tells you a direction; gradient descent is the rule that turns that direction into an update:

where θ stands for whichever parameter is being updated, and α (the learning rate) controls how big a step to take in that downhill direction. In code, for each parameter:

param -= alpha * grad

The choice of * here instead of @ matters and is easy to get backwards:

• * — element-wise (or scalar) multiplication. Each entry of grad is scaled by alpha independently, and the result has the same shape as param — exactly what an update needs.
• @ — matrix multiplication. It combines entries across dimensions and generally changes the shape of the result — not what you want when updating a parameter vector element-by-element.

Both parameters follow the identical rule — a weight vector update and a scalar bias update are the same formula, just with different shapes:

w -= alpha * dw
b -= alpha * db

Repeating "compute z → p → loss → gradient → update" for many iterations is the entire training loop.

12 Mistakes I Made

Click each one open. These are kept verbatim-in-spirit because the wrong version is often more instructive than the fix.

50000 <= income <= 10000          # impossible: lower bound > upper bound

income = np.random.randint(0, 10000, size=n)   # but the house-buying
                                                # rule assumed incomes
                                                # like 20,000–50,000

bias = np.random.randint(1, 10, size=10)

loss = -(y * p + (1 - y) * np.log(1 - p))   # missing np.log(p)

loss = -(
    y * np.log(p)
    + (1 - y) * np.log(1 - p)
)

p = np.clip(p, 1e-15, 1 - 1e-15)

self.data[:j]
error[j]

gradient /= len(param)   # should divide by number of examples

gradient /= len(data)

for j in range(len(param)):
    dw = gradient_for_feature_j   # overwritten every iteration

dw = []
for j in range(len(param)):
    gradient = ...
    dw.append(gradient)

param -= alpha @ grad

param -= alpha * grad

13 The Final Working Pipeline

Put together, with clipping for numerical safety and a loss history for later verification:

import numpy as np



class LogisticRegresssion:

    def __init__(self,data,output,param,bias):
        self.data = data
        self.output = output
        self.param = param
        self.bias = bias

    def calculate_z(self):
        z = self.data @ self.param + self.bias
        return z 

    def sigmoid(self):
        z = self.calculate_z()

        p = 1 / (1 + np.exp(-z))
        return p

    def calculate_loss(self):
        p = self.sigmoid()

        loss = - (self.output * np.log(p) + (1 - self.output) * np.log(1-p))

        return np.mean(loss)

    def calculate_gradient(self):
        p = self.sigmoid()
        y = self.output

        error = p - y
        dw = []
        for j in range(len(self.param)):
            grad = np.sum(error * self.data[:,j]) / len(self.data)
            dw.append(grad)

        db = np.sum(error)/len(self.data)
        return np.append(dw,db)


    def gradient_descent(self,alpha=0.1,iteration=1000):
        
        for i in range(iteration):
            loss = self.calculate_loss()
            print("loss_before: ",loss,end=" ")
            grad_val = self.calculate_gradient()

            self.param -= alpha * grad_val[:-1]
            self.bias -= alpha * grad_val[-1]

            loss = self.calculate_loss()
            print("param: ",self.param , " " , self.bias, "after_loss: ",loss)

        return self.param

14 How to Verify It Actually Works

1. Loss should trend down, then flatten

Under full-batch gradient descent with a sane learning rate, loss_history should look roughly like this — decreasing, with diminishing steps as it approaches a floor:

2. Inspect predictions, not just loss

p = model.predict_proba(X)
prediction = (p >= 0.5).astype(int)

3. Compute accuracy

accuracy = np.mean(prediction == y)

Loss and accuracy check different things, and relying on loss alone can be misleading. Loss keeps decreasing whenever probabilities move even slightly closer to the true labels — it's a continuous, fine-grained signal. Accuracy only cares whether the final thresholded class (0.5 cutoff) is right or wrong. A model can improve its loss for a while without its accuracy changing at all (probabilities shifting from 0.6 to 0.7 doesn't flip any predictions), so both checks belong in the verification routine, not just one.

15 What I Actually Learned

1. Logistic regression is not simply "linear regression + sigmoid" bolted on — the sigmoid changes what the output means, and that changes the loss function and the gradient too.
2. Sigmoid's only job is converting an unbounded linear score into something interpretable as a probability.
3. Log loss isn't an arbitrary formula — it falls out of maximum likelihood under a Bernoulli model of the labels.
4. The dataset supplies X and y; the model supplies (and learns) w and b, and computes z and p from them.
5. Loss measures how bad the current predictions are, right now.
6. Gradient tells you which direction reduces that loss — it's a slope, not a score.
7. The gradient is derived from the loss via calculus; the two are connected mathematically, but the numeric loss value is never plugged directly into the gradient computation.
8. Rows are training examples, columns are features — always, and every indexing bug in this project traced back to forgetting that.
9. m is the number of examples; n is the number of features. They are different numbers and get used in different places.
10. The parameter update is always \(\theta \leftarrow \theta - \alpha \nabla J\) — element-wise, using *, never @.
11. A decreasing loss curve is the first and most important sanity check that gradient descent is actually working.
12. Numerical symptoms like log(0) are usually messages about a deeper problem — bad feature scaling or an unstable learning rate — not just a clipping exercise.