Linear and logistic regression from scratch

Deriving the normal equations and gradient updates, implementing both models in NumPy, and verifying every implementation against scikit-learn.

Linear regression and the normal equations

Minimising ||Xw - y||^2 has a closed-form solution: w = (X^T X)^-1 X^T y. It is exact and instantaneous for small problems, and it is the reference you should check a gradient-descent implementation against.

import numpy as np

rng = np.random.default_rng(0)
n, d = 500, 4

X = rng.normal(size=(n, d))
X = np.c_[np.ones(n), X]                 # add an intercept column
true_w = np.array([1.5, -2.0, 0.5, 3.0, 0.0])
y = X @ true_w + 0.3 * rng.normal(size=n)

def normal_equations(X, y, ridge=0.0):
    A = X.T @ X + ridge * np.eye(X.shape[1])
    return np.linalg.solve(A, X.T @ y)   # solve, never invert

w_closed = normal_equations(X, y)
print(np.round(w_closed, 4))
print("MSE:", np.mean((X @ w_closed - y) ** 2))

# the same answer by gradient descent, which generalises to models with no closed form
w = np.zeros(d + 1)
for step in range(2000):
    grad = 2 * X.T @ (X @ w - y) / n
    w -= 0.1 * grad
np.allclose(w, w_closed, atol=1e-6)

Logistic regression and its gradient

With p = sigmoid(Xw) and binary cross-entropy, the derivative collapses to the same clean form as linear regression: X^T (p - y) / n. That is not a coincidence — it is what an exponential-family likelihood with a matching link gives you.

def sigmoid(z):
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))   # clip for stability

def fit_logistic(X, y, lr=0.5, steps=5000, l2=0.0):
    w = np.zeros(X.shape[1])
    losses = []
    for _ in range(steps):
        p = sigmoid(X @ w)
        eps = 1e-12
        loss = -(y * np.log(p + eps) + (1 - y) * np.log(1 - p + eps)).mean() + 0.5 * l2 * (w[1:] ** 2).sum()
        grad = X.T @ (p - y) / len(y) + l2 * np.r_[0.0, w[1:]]
        w -= lr * grad
        losses.append(loss)
    return w, losses

Xb = np.c_[np.ones(n), rng.normal(size=(n, 2))]
yb = (Xb @ np.array([0.5, 1.5, -1.0]) + 0.5 * rng.normal(size=n) > 0).astype(float)

w_log, losses = fit_logistic(Xb, yb)
p = sigmoid(Xb @ w_log)
acc = ((p > 0.5) == yb).mean()

from sklearn.linear_model import LogisticRegression
sk = LogisticRegression(penalty=None, max_iter=5000).fit(Xb[:, 1:], yb)
print(round(float(acc), 4))
print(np.round(w_log[1:], 3), np.round(sk.coef_[0], 3))

• Feature scaling matters much more here than in linear regression: uncentred features make the loss surface a narrow valley and the convergence painfully slow.
• The p - y gradient is already bounded in [-1, 1] per example, so logistic regression has no exploding-gradient problem — but it does have perfect separability, where weights grow without bound.
• Regularisation fixes separability and is essentially always appropriate when d is comparable to n.
• Accuracy hides calibration. Always inspect predicted probabilities, or a reliability curve, not just the thresholded decisions.

Verifying against scikit-learn

from sklearn.linear_model import LinearRegression, Ridge
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

# linear: compare coefficients directly
lr = LinearRegression().fit(X[:, 1:], y)
np.allclose(w_closed, np.r_[lr.intercept_, lr.coef_], atol=1e-8)

# ridge: sklearn scales the penalty by n/2, so match it exactly
alpha = 2.0
w_ridge = normal_equations(X, y, ridge=alpha * n / 2)
sk_ridge = Ridge(alpha=alpha, fit_intercept=True).fit(X[:, 1:], y)
np.allclose(w_ridge, np.r_[sk_ridge.intercept_, sk_ridge.coef_], atol=1e-6)

# logistic: sklearn regularises by default and applies C = 1 / lambda
Xs = StandardScaler().fit_transform(Xb[:, 1:])
Xstd = np.c_[np.ones(n), Xs]
w_std, _ = fit_logistic(Xstd, yb, lr=1.0, steps=8000, l2=1 / (2 * n))
sk2 = LogisticRegression(C=2 * n, max_iter=10000).fit(Xs, yb)
np.allclose(w_std[1:], sk2.coef_[0], atol=1e-2)

• Check the loss decreases monotonically before checking anything else; if it does not, the learning rate is too high or the gradient sign is wrong.
• Gradient-check one parameter group against a finite difference before trusting a from-scratch implementation.
• Match regularisation conventions before comparing coefficients: libraries differ by factors of n and 2, which is the most common source of a spurious mismatch.
• Evaluate on a held-out split. Fitting a synthetic dataset perfectly tells you the optimiser works, not that the model generalises.

FAQ

Related

Optimisation and gradient descent variants Information theory for machine learning

Last refreshed 2026-09-18.