Linear and logistic regression from scratch

Lesson 01 described a model as a program with blanks and training as a search for the numbers that fill them. This lesson does the search by hand for the two models that every other model is measured against. Linear regression predicts a number and has an exact formula. Logistic regression predicts a probability and does not, so it introduces gradient descent, the search that lesson 04 generalises and that every neural network uses. Each implementation is a few lines of numpy and is checked against scikit-learn, so you can see that the library is doing the same arithmetic.

Linear regression: the model and the loss

The model is a weighted sum. With one feature x, the prediction is w0 + w1 * x: an intercept and a slope. With d features it is the intercept plus one weight per feature. The parameters are the weights, and training chooses them.

The loss is the mean squared error, the average over examples of the squared difference between prediction and label. Squaring makes every error positive, punishes large errors much more than small ones, and gives a loss whose minimum can be found with algebra rather than search.

Work a tiny case in your head. Three examples with x equal to 1, 2 and 3 and labels 2, 3 and 5. The mean of x is 2 and the mean of the labels is 10 divided by 3, about 3.333. For a single feature, the best slope is the sum of (x minus its mean) times (label minus its mean), divided by the sum of (x minus its mean) squared. The numerator is (minus 1 times minus 1.333) plus 0 plus (1 times 1.667), which is 3.0. The denominator is 1 plus 0 plus 1, which is 2. The slope is 1.5. The intercept is the label mean minus slope times the x mean, 3.333 minus 3.0, which is 0.333. The line predicts 1.833, 3.333 and 4.833, missing the labels by 0.167, 0.333 and 0.167. No other line does better on the squared loss. That formula is the one-feature case of the general solution below.

The closed form

Stack the examples into a matrix X with one row per example and one column per feature, plus a column of ones so the intercept is handled like any other weight. Stack the labels into a vector y. The squared loss is then a bowl-shaped function of the weight vector w, and the bottom of a bowl is where its slope is zero in every direction. Setting the gradient to zero gives the normal equations: X transposed times X, times w, equals X transposed times y. That is a small square linear system, and np.linalg.solve finds w directly.

import numpy as np

rng = np.random.default_rng(0)
n = 50
x = rng.uniform(0, 10, size=n)              # hours of study, say
noise = rng.normal(0, 1.5, size=n)
y = 3.0 + 2.0 * x + noise                   # true line: intercept 3, slope 2

# Design matrix: a column of ones for the intercept, then x.
X = np.column_stack([np.ones(n), x])

# Normal equations: w = (X^T X)^-1 X^T y
w = np.linalg.solve(X.T @ X, X.T @ y)
print("closed-form intercept, slope:", np.round(w, 3))

residuals = y - X @ w
print("mean squared error:", round(float(np.mean(residuals ** 2)), 3))

closed-form intercept, slope: [2.636 2.079]
mean squared error: 2.235

The data was generated from intercept 3 and slope 2 with noise of standard deviation 1.5, so the fitted 2.636 and 2.079 are close but not exact. Fifty noisy points cannot pin down the truth, and the mean squared error of 2.235 is about the noise variance (1.5 squared is 2.25). That is what a well-fitted model on noisy data looks like: the residual error is the noise, and nothing more can be removed.

np.linalg.solve rather than computing an explicit inverse is a habit worth forming. Solving the system is faster and numerically safer than inverting the matrix and multiplying. For very wide feature sets or features that are near copies of each other, the matrix is close to singular and even solve struggles. Libraries use a decomposition instead, and add a small penalty on the weights, which is what “ridge regression” means. The classical ML hub covers that family.

Check against scikit-learn

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 50
x = rng.uniform(0, 10, size=n)
y = 3.0 + 2.0 * x + rng.normal(0, 1.5, size=n)
X = np.column_stack([np.ones(n), x])
w = np.linalg.solve(X.T @ X, X.T @ y)

sk = LinearRegression().fit(x.reshape(-1, 1), y)
print("numpy  :", np.round(w, 6))
print("sklearn:", np.round([sk.intercept_, sk.coef_[0]], 6))
print("max difference:", float(np.max(np.abs(w - [sk.intercept_, sk.coef_[0]]))))

numpy  : [2.636207 2.078504]
sklearn: [2.636207 2.078504]
max difference: 1.9539925233402755e-14

The two agree to fourteen decimal places, which is floating-point rounding. LinearRegression adds the intercept column for you, which is why it is given x as a column rather than X.

Classification: why squared error breaks

Now the label is yes or no, coded 1 or 0, as in the churn table of lesson 02. You could fit a line to the 0s and 1s and call anything above 0.5 a yes. It works badly for two reasons. The line’s output is unbounded, so for an extreme input it predicts 3.7 or minus 2, numbers that mean nothing as answers. And squared error punishes a confident correct answer: if the label is 1 and the line outputs 4, the loss is 9, so training pulls the line down towards examples it has already got right.

The fix has two parts. Keep the weighted sum, because it is the cheapest way to combine features, but push its output through a squashing function that maps any real number to the interval between 0 and 1. Then read the result as the probability that the label is 1, and choose a loss that rewards probabilities for being right rather than for being near a target number.

Logistic regression. A weighted sum of the features gives a score z, the sigmoid squashes z into a probability, and a threshold on the probability gives the class. Training adjusts the weights so that the probabilities match the labels.

The sigmoid

The squashing function is the sigmoid, 1 divided by (1 plus e to the minus z). At z equal to 0 it gives 0.5. At z equal to 2 it gives about 0.88, at 4 about 0.98, and at minus 2 about 0.12. Large positive scores become probabilities near 1, large negative scores near 0, and the transition is smooth, which matters because the search needs a slope to follow. The hidden rule that generated the churn labels in lesson 02 was exactly this: a weighted sum pushed through a sigmoid.

The log loss

For one example with label 1 and predicted probability p, the loss is minus the natural log of p. If p is 0.9 the loss is about 0.105. If p is 0.5 it is 0.693. If p is 0.1 it is 2.303. The loss rises slowly for small mistakes and steeply for confident wrong answers. For label 0 the same thing applies to (1 minus p). The mean over all examples is the log loss, also called binary cross-entropy. It is the loss that makes the predicted p behave like a calibrated probability, and it is the loss behind every language model’s next-token prediction, scaled up to a vocabulary of many thousands of choices.

The gradient

With the sigmoid and the log loss, something tidy happens. The derivative of the mean loss with respect to the weight vector is X transposed times (p minus y), divided by n. Read it as: for each example, the error is the predicted probability minus the true label, a number between minus 1 and 1, and each weight’s gradient is the average of that error times the corresponding feature value. Weights attached to features that are large when the model is wrong get the biggest nudge. There is no closed form because p depends on w through the sigmoid, so the equation “gradient equals zero” cannot be solved algebraically. Instead, step repeatedly in the direction opposite to the gradient.

import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def log_loss(y, p):
    eps = 1e-12
    p = np.clip(p, eps, 1 - eps)
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

rng = np.random.default_rng(0)
n = 200
# Two features. Class 1 tends to have larger x1 + x2.
X_raw = rng.normal(0, 1, size=(n, 2))
true_w = np.array([-0.5, 2.0, -1.5])          # bias, w1, w2
logits = true_w[0] + X_raw @ true_w[1:]
y = (rng.random(n) < sigmoid(logits)).astype(float)

X = np.column_stack([np.ones(n), X_raw])      # add the bias column
w = np.zeros(3)
lr = 0.5
for step in range(1, 501):
    p = sigmoid(X @ w)
    grad = X.T @ (p - y) / n                  # dL/dw for mean log loss
    w -= lr * grad
    if step in (1, 10, 50, 100, 500):
        print(f"step {step:4d}  loss {log_loss(y, p):.4f}  w {np.round(w, 3)}")

acc = float(np.mean((sigmoid(X @ w) >= 0.5) == y))
print("training accuracy:", round(acc, 3))

step    1  loss 0.6931  w [-0.04   0.134 -0.086]
step   10  loss 0.4789  w [-0.21   0.851 -0.555]
step   50  loss 0.3940  w [-0.369  1.765 -1.227]
step  100  loss 0.3861  w [-0.435  2.089 -1.485]
step  500  loss 0.3846  w [-0.488  2.319 -1.67 ]
training accuracy: 0.84

The first loss is 0.6931, which is the natural log of 2: with all weights at zero the model says 0.5 for everything, and that is what log loss costs when it has no opinion. By step 10 the weights have the right signs and the loss has dropped by a third. By step 100 they are close to the true values of minus 0.5, 2.0 and minus 1.5. Between steps 100 and 500 the loss barely moves while the weights keep growing a little, because with finite noisy data the fitted weights overshoot the generating ones slightly. Accuracy of 0.84 is near the ceiling for this data: the labels were drawn from probabilities, so some examples are on the wrong side of their own probability and no weights can recover them.

The clipping inside log_loss prevents taking the log of exactly 0, which numpy would report as negative infinity. Libraries do the same thing quietly.

Check against scikit-learn

import numpy as np
from sklearn.linear_model import LogisticRegression

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

rng = np.random.default_rng(0)
n = 200
X_raw = rng.normal(0, 1, size=(n, 2))
true_w = np.array([-0.5, 2.0, -1.5])
y = (rng.random(n) < sigmoid(true_w[0] + X_raw @ true_w[1:])).astype(float)

X = np.column_stack([np.ones(n), X_raw])
w = np.zeros(3)
for _ in range(5000):
    w -= 0.5 * X.T @ (sigmoid(X @ w) - y) / n

# C is the inverse regularisation strength. A huge C means "almost none",
# which is what our plain gradient descent is doing.
sk = LogisticRegression(C=1e6, max_iter=5000).fit(X_raw, y)
print("numpy  :", np.round(w, 3))
print("sklearn:", np.round(np.r_[sk.intercept_, sk.coef_[0]], 3))
print("true   :", true_w)
print("prediction agreement:",
      round(float(np.mean((sigmoid(X @ w) >= 0.5) == (sk.predict(X_raw) == 1))), 3))

numpy  : [-0.488  2.32  -1.671]
sklearn: [-0.488  2.32  -1.67 ]
true   : [-0.5  2.  -1.5]
prediction agreement: 1.0

The weights match to three decimals and every prediction agrees. Two things had to be set to get there. First, the numpy loop ran for 5,000 steps so that it had fully converged, since at 500 steps it was still creeping. Second, scikit-learn’s LogisticRegression applies an L2 penalty on the weights by default, controlled by C, which pulls them towards zero. The default C=1.0 would have given smaller weights than the plain loop. Setting C=1e6 switches the penalty off in effect. That default is sensible for real work, because it stops the weights running away on separable data, but it means the library is minimising a slightly different loss from the one on this page, and you should know that before comparing numbers.

Where linear models sit

Linear and logistic regression remain the first model to try on any tabular problem. They train in milliseconds, their weights can be read (a weight of 2.3 on support_calls means each extra call raises the churn score by 2.3), and they set a baseline that a fancier model must beat by a worthwhile margin. Their limit is the straight line: they cannot represent an interaction such as “high spend matters only for new customers” unless you build that feature by hand. Lesson 05 introduces trees, which discover interactions on their own, and lesson 07 stacks sigmoid-like units into a network that learns its own features. Both still end in the loss and gradient you have now seen.

Exercise

Add a second feature to the linear regression data, x2 = rng.uniform(0, 5, size=n), with true weight minus 1.0, and extend X to three columns. Confirm that np.linalg.solve recovers all three weights and that LinearRegression agrees.
In the logistic regression loop, try learning rates of 0.05 and 5.0 and print the loss at the same steps. Describe what each does, then read lesson 04 for the explanation.
Fit LogisticRegression with the default C=1.0 and compare its weights with the C=1e6 ones. Which direction does the penalty move them, and does the prediction agreement change.

Summary

Linear regression is a weighted sum with a squared-error loss. The loss is a bowl, and the normal equations find its bottom in one np.linalg.solve call.
For a yes-or-no label, the sigmoid turns the weighted sum into a probability and the log loss scores that probability, punishing confident mistakes steeply.
The gradient of the log loss is X transposed times (prediction minus label), averaged. There is no closed form, so gradient descent steps against that gradient until the loss stops falling.
The from-scratch versions match scikit-learn once regularisation is switched off and the loop has converged. The library adds a penalty, a faster solver and multiclass support, not a different model.

Linear and logistic regression from scratch

You will learn

Linear regression: the model and the loss

The closed form

Check against scikit-learn

Classification: why squared error breaks

The sigmoid

The log loss

The gradient

Check against scikit-learn

Where linear models sit

Exercise

Summary

Further reading