Logistic regression is used to model a binary label (i.e. 0/1) through a linear function followed by a sigmoid. Broadly speaking, the sigmoid is the common approach for modelling a Bernoulli distribution, which is a statistical model for predicting outcomes involving two distinct categories. Training a model involves maximizing conditional likelihood, or equivalently minimizing binary cross-entropy.

class LogisticRegression(ABC):
    """Base wrapper for the training loop."""

    def _sigmoid(self, logits: Tensor) -> Tensor:
        """Used to obtain probabilities from specified logits"""
        return 1/(1+torch.exp(-logits))

    def fit(self, X: Tensor, y: Tensor) -> LogisticRegression:
        for _ in range(self.n_iters):
            logits = X @ self.weights

            # Apply the sigmoid and the log likelihood function
            p = self._sigmoid(logits)
            p = torch.clamp(p, min=1e-12, max=1-1e-12)
            loss = -torch.mean(y*torch.log(p) + (1-y)*torch.log(1-p))

            # Update step
            step = self.optimizer(X, y)   # N
            self.weights = self.weights - self.lr * step

class GradientDescent(LogisticRegression):
    """Calculate gradient using gradient descent (first order method)."""

    def optimizer(self, X: Tensor, y: Tensor) -> Tensor:
        """Calculate the gradient"""
        logits = X @ self.weights   # (MxN) x (N) -> M
        p = self._sigmoid(logits)   # M
        return torch.mean(X.T * (p-y), dim=1) # (NxM) x (M) -> N x M -> N


class NewtonsMethod(LogisticRegression):
    """Calculate gradient using Newton's method (second order method)."""

    def optimizer(self, X: Tensor, y: Tensor) -> Tensor:
        """Calculate the gradient."""
        z = X @ self.weights    # (MxN) x N -> M
        p = self._sigmoid(z)    # M
        S = torch.diag(p*(1-p)) # MxM

        # Dimensions: (NxM x MxM x MxN)^-1 x NxM x M -> N
        H = X.T @ S @ X
        g = X.T @ (p - y)

        # Apply a dampening term so the system becomes stable
        eps = 1e-4
        I = torch.eye(H.shape[0], dtype=H.dtype, device=H.device)
        step = torch.linalg.solve(H + eps * I, g)
        return step

Model Setup

Setup. Let $X \in \mathbb{R}^{N \times p}$ contain rows $x_i^\top \in \mathbb{R}^{1 \times p}$ and let $y \in \{0,1\}^N$ contain the binary targets. If a bias is used, either augment each feature vector with a leading 1 and absorb the bias into $w$, or keep the bias separate. For coding convenience the former approach is preferred since it streamlines calculations (see here). The explanation below keeps the bias separate.

To obtain binary outputs, a common approach is to apply thresholding at $0.5$: $$ \hat y_i = \begin{cases} 1, & p_i > 0.5, \\ 0, & \text{otherwise}. \end{cases} $$

Objective

Optimization

Generally, first-order methods are commonly used for logistic regression. For smaller, more difficult problems, second-order methods such as Newton’s method and IRLS are also common because they can converge in fewer iterations, at the cost of being more compute-expensive.

First-Order Methods

Note that $X^\top(p-y) \in \mathbb{R}^p$, so the gradient matches the weights vector. If an $L_2$ penalty $\frac{\lambda}{2N}\lVert w\rVert_2^2$ is added, then $$ \nabla_w L_{\mathrm{reg}}(w,b) = \frac{1}{N}X^\top(p-y) + \frac{\lambda}{N}w, $$ while the bias derivative is commonly left unregularized.

Second-Order Methods

For Newton’s method and IRLS (iteratively reweighted least squares), use the augmented matrix $\tilde X = [\mathbf 1; X]$, the augmented parameter $\tilde w = (b,w)$, and the diagonal matrix $$ S := \operatorname{diag}(p_i(1-p_i)). $$

Implementation notes. Main implementation notes include:

• For stability, the predicted probabilities are clamped to $[10^{-12}, 1-10^{-12}]$ before taking logarithms.
• In Newton’s method, a damping term term is used, replacing $H$ with $H+\varepsilon I$ for $\varepsilon=10^{-4}$, to stabilize the linear solver when the Hessian is ill-conditioned (i.e. when the columns are nearly linearly dependent).
• Experiments confirm that Newton’s method converged faster and used $\eta=1.0$, while gradient descent still required a relatively large learning rate, $\eta=0.1$.

Citation

@misc{vasile2026logistic-regression-numpy,
    title = "Logistic Regression",
    author = "Razvan Florian Vasile",
    note = "Personal blog",
    year = "2026",
    month = "March",
    url = "https://rfvasile.github.io/posts/experiments/logistic_regression/"
}

References