Logistic Regression From Scratch in NumPy, With L2

Train p = σ(Xw + b) by batch gradient descent on

L = −(1/n) Σ [y_i log p_i + (1 − y_i) log(1 − p_i)] + (lam/2)‖w‖²,

starting from w = 0, b = 0, for steps iterations with learning rate lr. The bias is not penalised. Compute both gradients from the current parameters, then update both. Return [w₁, …, w_d, b]. Logits in the tests can exceed 700 in magnitude.

Examples

Input X = [[1], [-1]], y = [1, 0], lr = 1, steps = 1, lam = 0

Output [0.5, 0]

Both predictions start at ½. One step moves w to 0.5 and leaves b at 0.

Input X = [[0, 0], [1, 1], [2, 2], [3, 3]], y = [0, 0, 1, 1], lr = 0.5, steps = 500, lam = 0.1

Output [1.0112, 1.0112, -3.0334]

Separable data. The penalty keeps the weights finite.

Input X = [[1000], [-1000]], y = [1, 0], lr = 0.002, steps = 3, lam = 0

Output [1, 0]

After one step the logits are ±1000. A naive sigmoid overflows here.

Constraints

1 ≤ n ≤ 1000, 1 ≤ d ≤ 20
0 ≤ steps ≤ 2000
lam ≥ 0
Checked to a relative tolerance of 10⁻⁷

Hints

Hint 1

The derivative of the per-example loss with respect to the logit zᵢ is pᵢ − yᵢ. Start from there.

Hint 2

∂L/∂w = (1/n)Xᵀ(p − y) + lam·w and ∂L/∂b = mean(p − y). The penalty does not touch b.

Hint 3

For the sigmoid, compute 1/(1 + e^(−z)) where z ≥ 0 and e^z/(1 + e^z) where z < 0, so the exponential never sees a large positive number.

Where candidates go wrong

Penalising the bias. The tests do not, and it biases every prediction towards one half.
A sigmoid written as 1/(1 + np.exp(−z)): with z = −1000 the exponential overflows to inf. The answer, 0, is right, but NumPy warns and some setups treat warnings as errors.
Writing the gradient as Xᵀ(y − p) and then adding it. One sign error and the loss climbs.
Forgetting the 1/n on the data term but keeping lam·w unscaled, or the other way round.

The idea

For one example, ℓ = −y log σ(z) − (1 − y) log(1 − σ(z)) with z = x·w + b. Using σ′(z) = σ(z)(1 − σ(z)), the derivative collapses to ∂ℓ/∂z = σ(z) − y. The chain rule then gives ∂ℓ/∂w = (p − y)x and ∂ℓ/∂b = p − y. Averaging over the batch and adding the penalty's gradient lam·w:

g_w = Xᵀ(p − y)/n + lam·w, g_b = mean(p − y).

The loss is convex, so gradient descent with a small enough rate finds the unique minimum when lam > 0. Without the penalty and with separable data there is no minimum, and ‖w‖ grows without bound, which is the second example.

The sigmoid is computed on two branches so that exp is only ever called on non-positive numbers.

Time O(steps · nd) · Space O(nd)

Reference solution

import numpy as np

def sigmoid(z):
    # Two branches, so that exp only ever sees a non-positive argument
    out = np.empty_like(z)
    pos = z >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-z[pos]))
    ez = np.exp(z[~pos])
    out[~pos] = ez / (1.0 + ez)
    return out

def fit_logistic(X, y, lr, steps, lam):
    X = np.asarray(X, dtype=float)
    y = np.asarray(y, dtype=float)
    n, d = X.shape
    w, b = np.zeros(d), 0.0
    for _ in range(steps):
        err = sigmoid(X @ w + b) - y        # dL/dz for each example (before the 1/n)
        gw = X.T @ err / n + lam * w        # the penalty does not touch the bias
        gb = err.mean()
        w -= lr * gw
        b -= lr * gb
    return [*w.tolist(), float(b)]

Console⌘/Ctrl + Enter runs

Python runs in your browser with NumPy. The first run downloads the runtime, which takes a few seconds.

Why interviewers ask it

This is the most common "implement a model" question in ML interviews. It tests the gradient of cross-entropy, whether the bias is penalised, and whether your sigmoid survives a logit of −1000. The derivation is usually asked before any code.

The theory behind it is in the lesson Linear and logistic regression from scratch.

Keep going

Easy10 min

Logistic regression with an L2 penalty