diff --git a/README.md b/README.md
index beec699..1675380 100644
--- a/README.md
+++ b/README.md
@@ -21,11 +21,16 @@ Every deep learning framework is, at its core, a graph that records operations a
 - Minibatch stochastic gradient descent with feature standardization for conditioning
 - The equivalence of the closed-form and iterative solutions on a convex loss
 - Pseudoinverse fallback for rank-deficient design matrices
+- Binary logistic regression as a Bernoulli likelihood, trained by minimizing cross-entropy
+- The `sigmoid(z) - y` gradient identity, the same residual form that least squares has
+- Numerically stable sigmoid (sign-branched) and log loss via `softplus(z) - y·z` (`logaddexp`)
+- The linear decision boundary as the level set where the logit is zero
 
 ## What's implemented
 
 - **Gradient descent + autodiff-lite**: a scalar reverse-mode autograd engine (`Value`) with `+`, `*`, `**`, division, and `relu`/`tanh`/`exp`/`log`/`sigmoid`, plus an `SGD` optimizer (with momentum) and a `minimize` training loop. Gradients are checked against finite differences and against numpy for softmax cross-entropy.
 - **Linear regression, two ways**: `fit_normal_equation` solves ordinary least squares exactly by setting the gradient to zero and solving `(AᵀA + λR)θ = Aᵀy`, with optional ridge and a least-norm fallback when `AᵀA` is singular. `fit_sgd` fits the same model with minibatch gradient descent on standardized features, then maps the weights back to raw feature space. Both are checked to recover the true coefficients and to agree with each other, so the exact-vs-iterative tradeoff is measurable.
+- **Logistic regression + cross-entropy, decision boundary demo**: `src/logreg.py` trains a binary classifier by minibatch SGD, using the fact that the gradient of the mean cross-entropy with respect to the logits is exactly `sigmoid(z) - y`, the same residual form linear regression has. The sigmoid branches on the sign of the logit so `exp` never overflows, and the loss is computed as `softplus(z) - y·z` through `numpy.logaddexp` so a confidently-wrong prediction gives a large finite loss instead of `inf`. `decision_boundary` returns the straight line where the model sits at 50 percent for a two-feature problem, the level set `w·x + b = 0`.
 
 ## Usage
 
@@ -55,6 +60,25 @@ print(exact.weights, exact.bias)   # ~[2, -3, 0.5], 1.5
 print(mse(approx, X, y))           # ~0, matches the closed form
 ```
 
+Fit a binary classifier and read its decision boundary:
+
+```python
+import numpy as np
+from src.logreg import fit_sgd, accuracy, decision_boundary
+
+rng = np.random.default_rng(0)
+X = rng.uniform(-3, 3, size=(400, 2))
+y = (X @ np.array([3.0, -2.0]) + 0.5 > 0).astype(int)
+
+model, history = fit_sgd(X, y, lr=0.2, epochs=300)
+print(accuracy(model, X, y))          # ~1.0 on separable data
+print(model.predict_proba(X[:3]))     # calibrated probabilities in (0, 1)
+
+# the 50 percent line, ready to plot against the two features
+xs = np.linspace(-3, 3, 50)
+boundary = decision_boundary(model, xs)
+```
+
 Differentiate an arbitrary scalar expression:
 
 ```python
diff --git a/src/logreg.py b/src/logreg.py
new file mode 100644
index 0000000..b7953fb
--- /dev/null
+++ b/src/logreg.py
@@ -0,0 +1,155 @@
+"""Binary logistic regression trained with cross-entropy loss.
+
+Logistic regression models P(y=1 | x) as sigmoid(w·x + b). Its loss is the
+negative log-likelihood of the Bernoulli labels, which is convex, so gradient
+descent finds the global optimum. The gradient has the same shape as linear
+regression: the residual (predicted probability minus label) times the input.
+The one thing you cannot do naively is the loss itself, because sigmoid and log
+overflow on large logits, so both the loss and the sigmoid are written in their
+numerically stable forms here. The decision boundary is the line where the
+logit is zero, and for two features that is a straight line you can plot.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import NDArray
+
+Array = NDArray[np.float64]
+
+
+@dataclass(frozen=True)
+class LogisticModel:
+    weights: Array
+    bias: float
+
+    def decision_function(self, X: Array) -> Array:
+        """Raw logits w·x + b, before the sigmoid."""
+        return _as_2d(X) @ self.weights + self.bias
+
+    def predict_proba(self, X: Array) -> Array:
+        return _sigmoid(self.decision_function(X))
+
+    def predict(self, X: Array, threshold: float = 0.5) -> Array:
+        return (self.predict_proba(X) >= threshold).astype(np.int64)
+
+
+def _sigmoid(z: Array) -> Array:
+    # branch on sign so neither exp overflows: exp(-|z|) is always in (0, 1]
+    out = np.empty_like(z, dtype=np.float64)
+    pos = z >= 0.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 bce_loss(model: LogisticModel, X: Array, y: Array) -> float:
+    """Mean binary cross-entropy in logits, via softplus(z) - y·z.
+
+    softplus(z) = log(1 + e^z) computed as logaddexp(0, z) never overflows, and
+    the identity BCE = softplus(z) - y·z avoids ever calling log on a sigmoid
+    that has saturated to exactly 0 or 1.
+    """
+    Xm, yv = _check_xy(X, y)
+    z = Xm @ model.weights + model.bias
+    return float(np.mean(np.logaddexp(0.0, z) - yv * z))
+
+
+def decision_boundary(model: LogisticModel, x1: Array) -> Array:
+    """For a 2-feature model, the x2 values where P(y=1)=0.5 along given x1.
+
+    The boundary is w0·x1 + w1·x2 + b = 0, so x2 = -(b + w0·x1) / w1. Only
+    defined when the model has exactly two features and w1 is nonzero (a
+    horizontal boundary has no functional form in x2).
+    """
+    if model.weights.shape != (2,):
+        raise ValueError("decision_boundary is only defined for two features")
+    w0, w1 = float(model.weights[0]), float(model.weights[1])
+    if w1 == 0.0:
+        raise ValueError("boundary is vertical in x2; w1 is zero")
+    x1v = np.asarray(x1, dtype=np.float64)
+    return -(model.bias + w0 * x1v) / w1
+
+
+def accuracy(model: LogisticModel, X: Array, y: Array) -> float:
+    Xm, yv = _check_xy(X, y)
+    return float(np.mean(model.predict(Xm) == yv))
+
+
+def _as_2d(X: Array) -> Array:
+    arr = np.asarray(X, dtype=np.float64)
+    return arr.reshape(-1, 1) if arr.ndim == 1 else arr
+
+
+def _check_xy(X: Array, y: Array) -> tuple[Array, Array]:
+    Xm = _as_2d(X)
+    yv = np.asarray(y, dtype=np.float64).reshape(-1)
+    if Xm.shape[0] == 0:
+        raise ValueError("need at least one sample")
+    if Xm.shape[0] != yv.shape[0]:
+        raise ValueError(f"X has {Xm.shape[0]} rows but y has {yv.shape[0]} entries")
+    if not np.all((yv == 0.0) | (yv == 1.0)):
+        raise ValueError("labels must be binary 0/1")
+    return Xm, yv
+
+
+def fit_sgd(
+    X: Array,
+    y: Array,
+    lr: float = 0.1,
+    epochs: int = 200,
+    batch_size: int = 32,
+    l2: float = 0.0,
+    seed: int = 0,
+) -> tuple[LogisticModel, list[float]]:
+    """Fit logistic regression by minibatch gradient descent on standardized X.
+
+    The gradient of the mean cross-entropy w.r.t. the logits is (p - y), the
+    same residual form as least squares, so a batch update is (1/m) Xᵀ(p - y).
+    Features are standardized to unit variance first for conditioning, then the
+    weights are mapped back to raw feature space so the returned model consumes
+    unscaled X. L2 shrinks the weights but not the bias. Returns the model and
+    the per-epoch training cross-entropy.
+    """
+    if lr <= 0.0:
+        raise ValueError("lr must be positive")
+    if epochs <= 0:
+        raise ValueError("epochs must be positive")
+    if batch_size <= 0:
+        raise ValueError("batch_size must be positive")
+    if l2 < 0.0:
+        raise ValueError("l2 must be non-negative")
+
+    Xm, yv = _check_xy(X, y)
+    n, d = Xm.shape
+
+    mu = Xm.mean(axis=0)
+    sigma = Xm.std(axis=0)
+    sigma = np.where(sigma == 0.0, 1.0, sigma)  # constant features carry no info
+    Xs = (Xm - mu) / sigma
+
+    rng = np.random.default_rng(seed)
+    w = np.zeros(d)
+    b = 0.0
+    history: list[float] = []
+
+    for _ in range(epochs):
+        order = rng.permutation(n)
+        for start in range(0, n, batch_size):
+            idx = order[start : start + batch_size]
+            xb, yb = Xs[idx], yv[idx]
+            m = xb.shape[0]
+            resid = _sigmoid(xb @ w + b) - yb
+            grad_w = (xb.T @ resid) / m + l2 * w
+            grad_b = float(resid.sum()) / m
+            w -= lr * grad_w
+            b -= lr * grad_b
+        z_full = Xs @ w + b
+        history.append(float(np.mean(np.logaddexp(0.0, z_full) - yv * z_full)))
+
+    weights = w / sigma
+    bias = b - float(np.sum(w * mu / sigma))
+    return LogisticModel(weights=weights, bias=bias), history
diff --git a/tests/test_logreg.py b/tests/test_logreg.py
new file mode 100644
index 0000000..6f5f6a5
--- /dev/null
+++ b/tests/test_logreg.py
@@ -0,0 +1,174 @@
+import dataclasses
+
+import numpy as np
+import pytest
+
+from src.logreg import (
+    LogisticModel,
+    accuracy,
+    bce_loss,
+    decision_boundary,
+    fit_sgd,
+)
+
+
+def make_blobs(n=400, w=(3.0, -2.0), b=0.5, seed=0):
+    """Two features, labels drawn from the true logistic model."""
+    rng = np.random.default_rng(seed)
+    X = rng.uniform(-3.0, 3.0, size=(n, len(w)))
+    z = X @ np.array(w) + b
+    p = 1.0 / (1.0 + np.exp(-z))
+    y = (rng.uniform(size=n) < p).astype(np.int64)
+    return X, y
+
+
+def separable_data(n=200, seed=1):
+    rng = np.random.default_rng(seed)
+    pos = rng.uniform(1.0, 3.0, size=(n // 2, 2))
+    neg = rng.uniform(-3.0, -1.0, size=(n // 2, 2))
+    X = np.vstack([pos, neg])
+    y = np.concatenate([np.ones(n // 2), np.zeros(n // 2)]).astype(np.int64)
+    return X, y
+
+
+def test_sigmoid_is_stable_on_extreme_logits():
+    model = LogisticModel(weights=np.array([1.0]), bias=0.0)
+    huge = np.array([-1000.0, 1000.0])
+    p = model.predict_proba(huge)
+    assert p[0] == pytest.approx(0.0, abs=1e-12)
+    assert p[1] == pytest.approx(1.0, abs=1e-12)
+    assert np.all(np.isfinite(p))
+
+
+def test_bce_matches_naive_formula_in_safe_range():
+    X, y = make_blobs(n=50, seed=3)
+    model = LogisticModel(weights=np.array([0.4, -0.6]), bias=0.1)
+    p = model.predict_proba(X)
+    naive = -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
+    assert bce_loss(model, X, y) == pytest.approx(naive, abs=1e-9)
+
+
+def test_bce_finite_on_saturated_logits():
+    # a confident-but-wrong prediction would make the naive log blow up
+    model = LogisticModel(weights=np.array([50.0]), bias=0.0)
+    X = np.array([[1.0]])
+    y = np.array([0])
+    loss = bce_loss(model, X, y)
+    assert np.isfinite(loss)
+    assert loss == pytest.approx(50.0, abs=1e-6)  # softplus(50) - 0 ≈ 50
+
+
+def test_recovers_true_weights_direction():
+    X, y = make_blobs(n=4000, w=(3.0, -2.0), b=0.5, seed=7)
+    model, history = fit_sgd(X, y, lr=0.2, epochs=400, batch_size=64)
+    # logistic weights are only identified up to the noise, so check direction
+    est = np.append(model.weights, model.bias)
+    true = np.array([3.0, -2.0, 0.5])
+    cos = est @ true / (np.linalg.norm(est) * np.linalg.norm(true))
+    assert cos > 0.98
+    assert history[-1] < history[0]
+
+
+def test_perfectly_separable_data_is_classified_perfectly():
+    X, y = separable_data(n=300, seed=2)
+    model, _ = fit_sgd(X, y, lr=0.3, epochs=300, batch_size=32)
+    assert accuracy(model, X, y) == pytest.approx(1.0)
+
+
+def test_loss_decreases_monotone_enough():
+    X, y = make_blobs(n=1000, seed=5)
+    _, history = fit_sgd(X, y, lr=0.2, epochs=200, batch_size=64)
+    assert len(history) == 200
+    assert history[-1] < 0.9 * history[0]
+
+
+def test_l2_shrinks_weights():
+    X, y = separable_data(n=300, seed=4)
+    plain, _ = fit_sgd(X, y, lr=0.2, epochs=200, l2=0.0)
+    ridged, _ = fit_sgd(X, y, lr=0.2, epochs=200, l2=2.0)
+    assert np.linalg.norm(ridged.weights) < np.linalg.norm(plain.weights)
+
+
+def test_predict_uses_threshold():
+    model = LogisticModel(weights=np.array([1.0]), bias=0.0)
+    X = np.array([[0.2]])  # sigmoid(0.2) ≈ 0.55
+    assert model.predict(X)[0] == 1
+    assert model.predict(X, threshold=0.6)[0] == 0
+
+
+def test_decision_boundary_line():
+    model = LogisticModel(weights=np.array([2.0, -4.0]), bias=1.0)
+    x1 = np.array([-1.0, 0.0, 1.0])
+    x2 = decision_boundary(model, x1)
+    # every (x1, x2) on the line must sit exactly at logit zero
+    pts = np.column_stack([x1, x2])
+    np.testing.assert_allclose(model.decision_function(pts), 0.0, atol=1e-9)
+    np.testing.assert_allclose(model.predict_proba(pts), 0.5, atol=1e-9)
+
+
+def test_decision_boundary_rejects_non_2d_models():
+    with pytest.raises(ValueError):
+        decision_boundary(LogisticModel(weights=np.array([1.0]), bias=0.0), [0.0])
+
+
+def test_decision_boundary_rejects_horizontal():
+    with pytest.raises(ValueError):
+        decision_boundary(LogisticModel(weights=np.array([1.0, 0.0]), bias=0.0), [0.0])
+
+
+def test_predict_proba_shape_and_range():
+    X, y = make_blobs(n=30, seed=6)
+    model, _ = fit_sgd(X, y, lr=0.1, epochs=50)
+    p = model.predict_proba(X)
+    assert p.shape == (30,)
+    assert np.all((p >= 0.0) & (p <= 1.0))
+
+
+def test_constant_feature_does_not_break_standardization():
+    X, y = separable_data(n=100, seed=8)
+    X = np.hstack([X, np.ones((X.shape[0], 1))])  # zero-variance column
+    model, _ = fit_sgd(X, y, lr=0.2, epochs=200)
+    assert np.all(np.isfinite(model.weights))
+    assert accuracy(model, X, y) > 0.95
+
+
+def test_single_feature_1d_input():
+    rng = np.random.default_rng(9)
+    x = rng.uniform(-3.0, 3.0, size=200)
+    y = (x > 0).astype(np.int64)
+    model, _ = fit_sgd(x, y, lr=0.3, epochs=300)
+    assert model.weights.shape == (1,)
+    assert model.predict(np.array([2.0]))[0] == 1
+    assert model.predict(np.array([-2.0]))[0] == 0
+
+
+def test_empty_input_raises():
+    with pytest.raises(ValueError):
+        fit_sgd(np.zeros((0, 2)), np.zeros(0))
+
+
+def test_mismatched_lengths_raise():
+    with pytest.raises(ValueError):
+        fit_sgd(np.zeros((5, 2)), np.zeros(4))
+
+
+def test_non_binary_labels_raise():
+    X = np.zeros((4, 2))
+    with pytest.raises(ValueError):
+        fit_sgd(X, np.array([0, 1, 2, 1]))
+    model = LogisticModel(weights=np.zeros(2), bias=0.0)
+    with pytest.raises(ValueError):
+        bce_loss(model, X, np.array([0, 1, 0, 5]))
+
+
+def test_invalid_hyperparams_raise():
+    X, y = make_blobs(n=20)
+    for kwargs in ({"lr": 0.0}, {"epochs": 0}, {"batch_size": 0}, {"l2": -1.0}):
+        with pytest.raises(ValueError):
+            fit_sgd(X, y, **kwargs)
+
+
+def test_logistic_model_is_frozen():
+    m = LogisticModel(weights=np.array([1.0]), bias=0.0)
+    with pytest.raises(dataclasses.FrozenInstanceError):
+        m.bias = 2.0