The essence of deep learning modelling

import numpy as np

def relu(z):    return np.maximum(0, z)
def sigmoid(z): return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

def forward(X, params):
    """X: [n_samples, n_features]. Cache everything the backward pass needs."""
    W1, b1, W2, b2 = params["W1"], params["b1"], params["W2"], params["b2"]

    Z1 = X @ W1 + b1          # linear:     [n, hidden]
    A1 = relu(Z1)             # non-linear: [n, hidden]
    Z2 = A1 @ W2 + b2         # linear:     [n, out]
    A2 = sigmoid(Z2)          # output:     [n, out]

    return A2, {"X": X, "Z1": Z1, "A1": A1, "Z2": Z2, "A2": A2}

rng = np.random.default_rng(0)
n_in, n_hidden, n_out = 2, 8, 1
params = {
    # He initialisation: scale by sqrt(2/fan_in) to keep the variance of
    # activations stable through ReLU, which discards half the distribution.
    "W1": rng.normal(0, np.sqrt(2 / n_in), (n_in, n_hidden)),
    "b1": np.zeros(n_hidden),
    "W2": rng.normal(0, np.sqrt(2 / n_hidden), (n_hidden, n_out)),
    "b2": np.zeros(n_out),
}

# XOR — not linearly separable, so a hidden layer is required.
X = np.array([[0., 0.], [0., 1.], [1., 0.], [1., 1.]])
y = np.array([[0.], [1.], [1.], [0.]])
predictions, cache = forward(X, params)

import numpy as np

def mse(y_pred, y_true):
    return float(np.mean((y_pred - y_true) ** 2))

def binary_cross_entropy(y_pred, y_true, eps=1e-12):
    p = np.clip(y_pred, eps, 1 - eps)          # clip, or log(0) = -inf
    return float(-np.mean(y_true * np.log(p) + (1 - y_true) * np.log(1 - p)))

y_true = np.array([[1.0]])
for p in [0.9, 0.5, 0.1, 0.01]:
    pred = np.array([[p]])
    print(f"p={p:<5} mse {mse(pred, y_true):.4f}   bce {binary_cross_entropy(pred, y_true):.4f}")
# MSE tops out at 1.0 no matter how wrong the prediction is.
# BCE grows without bound, which is the incentive you want.

# The gradient argument, which is the real reason:
#   MSE + sigmoid  ->  dL/dz = (a - y) * a * (1 - a)
#   BCE + sigmoid  ->  dL/dz = (a - y)
# The sigmoid derivative a(1-a) approaches 0 when a is near 0 or 1, so MSE's
# gradient vanishes precisely when the model is confidently wrong. BCE cancels
# that term exactly — which is why the pairing matters.

import numpy as np

def backward(y_true, cache, params):
    n = y_true.shape[0]
    X, A1, Z1, A2 = cache["X"], cache["A1"], cache["Z1"], cache["A2"]

    # Output layer. For sigmoid + binary cross-entropy this simplifies to
    # (A2 - y) — the sigmoid derivative cancels exactly.
    dZ2 = (A2 - y_true) / n                  # [n, out]
    dW2 = A1.T @ dZ2                         # [hidden, out]
    db2 = dZ2.sum(axis=0)

    # Propagate backwards through W2, then through the ReLU.
    dA1 = dZ2 @ params["W2"].T               # [n, hidden]
    dZ1 = dA1 * (Z1 > 0)                     # ReLU'(z) = 1 if z > 0 else 0
    dW1 = X.T @ dZ1
    db1 = dZ1.sum(axis=0)

    return {"W1": dW1, "b1": db1, "W2": dW2, "b2": db2}

# Gradient checking: compare against finite differences. Do this whenever you
# implement a gradient by hand — it catches every sign and transpose error.
def gradient_check(params, X, y, key, h=1e-6):
    analytic = backward(y, forward(X, params)[1], params)[key]
    numeric = np.zeros_like(params[key])
    it = np.nditer(params[key], flags=["multi_index"])
    while not it.finished:
        i = it.multi_index
        original = params[key][i]
        params[key][i] = original + h
        loss_plus = binary_cross_entropy(forward(X, params)[0], y)
        params[key][i] = original - h
        loss_minus = binary_cross_entropy(forward(X, params)[0], y)
        params[key][i] = original
        numeric[i] = (loss_plus - loss_minus) / (2 * h)
        it.iternext()
    rel = np.abs(analytic - numeric).max() / (np.abs(analytic).max() + 1e-12)
    return rel        # below ~1e-6 means the gradient is right

print("relative error:", gradient_check(params, X, y, "W1"))

# The full training loop, now that all three pieces exist.
for step in range(5000):
    preds, cache = forward(X, params)
    grads = backward(y, cache, params)
    for k in params:
        params[k] -= 0.5 * grads[k]
    if step % 1000 == 0:
        print(step, round(binary_cross_entropy(preds, y), 4))
print(forward(X, params)[0].round(3))     # solves XOR

import torch, torch.nn as nn

X = torch.tensor([[0., 0.], [0., 1.], [1., 0.], [1., 1.]])
y = torch.tensor([[0.], [1.], [1.], [0.]])

model = nn.Sequential(
    nn.Linear(2, 8),
    nn.ReLU(),
    nn.Linear(8, 1),
)
# BCEWithLogitsLoss takes RAW logits and applies the sigmoid internally, using
# the stable log-sum-exp form. Never add a Sigmoid layer and then use BCELoss.
loss_fn = nn.BCEWithLogitsLoss()
opt = torch.optim.Adam(model.parameters(), lr=0.05)

for step in range(2000):
    logits = model(X)                # 1. forward
    loss = loss_fn(logits, y)        # 2. cost

    opt.zero_grad()                  # 3a. clear old gradients — they ACCUMULATE
    loss.backward()                  # 3b. backward: autograd does the chain rule
    opt.step()                       # 3c. update

    if step % 500 == 0:
        print(step, round(loss.item(), 4))

with torch.no_grad():
    print(torch.sigmoid(model(X)).round())

# Inspecting what autograd computed — the same numbers the hand-written
# backward produced:
loss = loss_fn(model(X), y)
loss.backward()
for name, p in model.named_parameters():
    print(f"{name:16} grad norm {p.grad.norm().item():.4f}")