Fitting the Model

01 Learning objectives

02 Motivating example

In T-02 you were handed a logistic regression with coefficients already attached:

z = −3 + 1.2 × hours

and used it to predict exam pass probabilities. In T-03 you took another set of coefficients — z = −8 + 0.3x₁ + 0.5x₂ — and derived the decision boundary. In both cases the coefficients were given.

This tutorial asks the obvious follow-up question: where did those numbers come from?

For linear regression, from Unit 4, the answer was clean. The coefficients are the values that minimise the sum of squared residuals, and they can be written down in closed form as β = (XᵀX)⁻¹Xᵀy. Substitute the data, evaluate a matrix expression, and you have the answer. No iteration, no tuning, no ambiguity.

It is natural to hope that logistic regression works the same way. Fit a line by least squares, then squash it through the sigmoid — or perhaps minimise the squared error between the predicted probability and the binary label. Either approach seems plausible, and both are wrong.

Two things go wrong. First, least squares is not the right criterion for a binary outcome — it treats a confident mistake and a cautious mistake as equally bad in a way that makes no sense for probabilities. Second, even with the right criterion, there is no closed form. The equations that define the optimal coefficients cannot be solved algebraically; they must be solved numerically, by iteration.

This tutorial works through both issues. You will see what the right criterion is, why it is the right one, what its gradient looks like, why an algebraic solution is impossible, and how the numerical methods that do the job actually work.

03 Concept explanation

3.1 Why least squares is the wrong criterion

Suppose we try to fit logistic regression by minimising the sum of squared differences between the predicted probability and the binary label:

J(β) = Σ ( y_i − σ(z_i) )²

This is sometimes called the "squared error" or "Brier score" criterion, and it is not absurd — it is a legitimate way to measure the quality of probabilistic predictions. But it has two problems as a fitting criterion for logistic regression.

It penalises confident mistakes weakly. Consider two predictions for a case whose true label is 1:

• Model A predicts p = 0.9. Squared error: (1 − 0.9)² = 0.01.
• Model B predicts p = 0.6. Squared error: (1 − 0.6)² = 0.16.

Now consider two confident mistakes:

• Model C predicts p = 0.1. Squared error: (1 − 0.1)² = 0.81.
• Model D predicts p = 0.01. Squared error: (1 − 0.01)² = 0.98.

The difference between C and D is only 0.17, yet D is ten times more confident in a wrong answer. Squared error grows towards a ceiling of 1 and cannot express the fact that being confidently wrong is qualitatively worse than being cautiously wrong.

It is not the criterion the model assumes. Logistic regression models the probability that y = 1 as a Bernoulli random variable with parameter p. The natural way to fit a probabilistic model is by maximum likelihood, which asks: which parameter values make the observed data most probable? Minimising squared error is not the same as maximising likelihood, and the two produce different coefficient estimates.

3.2 Maximum likelihood

The principle of maximum likelihood is straightforward:

Choose the parameter values that make the observed data most probable under the assumed model.

Applied to logistic regression, this means finding the coefficients β that maximise the probability of the observed labels y₁, …, yₙ, given the observed features x₁, …, xₙ.

The likelihood function L(β) is exactly this probability, viewed as a function of the parameters. We then find the β that makes it as large as possible.

3.3 The likelihood for a single observation

For a single case with true label y ∈ {0, 1} and predicted probability p = σ(z), the probability of observing that label is:

• If y = 1: the probability is p.
• If y = 0: the probability is 1 − p.

These two cases can be written as a single expression:

P(y | x, β) = p^y · (1 − p)^(1−y)

Check the two cases. When y = 1, the expression becomes p¹ · (1 − p)⁰ = p. When y = 0, it becomes p⁰ · (1 − p)¹ = 1 − p. This compact form — sometimes called the Bernoulli likelihood — is worth recognising.

3.4 The log-likelihood and the log loss

Because the observations are assumed independent, the likelihood of the whole dataset is the product of the individual likelihoods:

L(β) = Π p_i^(y_i) · (1 − p_i)^(1−y_i)

Products are awkward to work with, and they underflow quickly on a computer. Taking the logarithm turns the product into a sum without changing the location of the maximum, because the logarithm is strictly increasing. This gives the log-likelihood:

ℓ(β) = Σ [ y_i log(p_i) + (1 − y_i) log(1 − p_i) ]

We maximise ℓ(β). Equivalently, we can minimise its negative:

J(β) = − Σ [ y_i log(p_i) + (1 − y_i) log(1 − p_i) ]

This quantity has several names: the negative log-likelihood, the log loss, the cross-entropy loss, or simply the logistic loss. They all refer to the same thing. Most machine learning libraries, including scikit-learn, minimise a version of this quantity (usually divided by n to give a mean rather than a sum).

The log loss behaves very differently from squared error. For a single observation, the contribution is:

• If y = 1: −log(p). This goes to +∞ as p → 0, and to 0 as p → 1.
• If y = 0: −log(1 − p). This goes to +∞ as p → 1, and to 0 as p → 0.

A confident mistake — y = 1 with p = 0.001 — incurs a loss of −log(0.001) ≈ 6.9, and grows without bound as the prediction becomes more confidently wrong. This is exactly the behaviour squared error could not express, and it is why the log loss is the right criterion for probabilistic classifiers.

3.5 The gradient

To find the maximum of ℓ(β), we need its gradient with respect to each coefficient. The derivation uses three facts:

1. The derivative of the sigmoid is σ′(z) = σ(z)(1 − σ(z)).
2. The chain rule applies through z_i = x_iᵀβ.
3. For a single observation, the log-likelihood simplifies dramatically.

Working through the calculus for a single observation:

∂ℓ_i/∂β_j = (y_i − p_i) · x_ij

Summing over all observations:

∂ℓ/∂β_j = Σ (y_i − p_i) · x_ij

In matrix form, with X the n × p design matrix, y the vector of labels, and p the vector of predicted probabilities:

∇ℓ(β) = Xᵀ (y − p)

This is a remarkably clean result, and it should look familiar. In linear regression, the gradient of the residual sum of squares is −2Xᵀ(y − Xβ). The structural form is identical: a matrix transpose times a vector of residuals. The only difference is that where linear regression uses the fitted values Xβ, logistic regression uses the fitted probabilities p = σ(Xβ).

3.6 Why there is no closed-form solution

Setting the gradient to zero gives the normal equations for logistic regression:

Xᵀ (y − p) = 0

where p = σ(Xβ) depends non-linearly on β. Expanding:

Xᵀ y = Xᵀ p = Xᵀ σ(Xβ)

The left-hand side is a linear function of the data. The right-hand side is a non-linear function of β, because σ is non-linear. There is no algebraic rearrangement that isolates β.

This is worth contrasting with ordinary least squares. In linear regression, the same gradient equation is Xᵀ(y − Xβ) = 0, which rearranges to XᵀXβ = Xᵀy — a linear system in β, solvable in one step. The difference is entirely due to the sigmoid: it introduces a non-linearity that cannot be undone algebraically.

The consequence is that logistic regression coefficients must be found by iterative numerical optimisation. This is not a limitation of the software; it is a mathematical fact about the model.

3.7 Gradient ascent

The simplest way to find the maximum is gradient ascent: start with an initial guess, then repeatedly step in the direction of the gradient, which points uphill.

β ← β + α · ∇ℓ(β)

where α > 0 is the learning rate (or step size). The algorithm terminates when the gradient is close enough to zero, or after a fixed number of iterations.

Three practical considerations:

• Choice of α. Too small and convergence is glacial; too large and the iteration can oscillate or diverge. There is no universal safe value; it depends on the data.
• Scaling. Features on different scales make the problem ill-conditioned, and gradient ascent becomes slow or unstable. Standardising features is almost always necessary in practice.
• Stopping rule. The iteration is usually stopped when the change in the log-likelihood (or the magnitude of the gradient) falls below a tolerance.

Gradient ascent is conceptually simple and easy to implement, but it can be slow. For logistic regression there is a much better method.

3.8 Newton–Raphson and IRLS

Newton's method uses second-order information — the curvature of the function — to take better steps. For maximisation, the update is:

β ← β − H⁻¹ ∇ℓ

where H is the Hessian matrix of second derivatives. For logistic regression, the Hessian has a clean form:

H = − Xᵀ W X

where W is the diagonal matrix with entries w_i = p_i(1 − p_i). Substituting into the update:

β ← β + (Xᵀ W X)⁻¹ Xᵀ (y − p)

This is the Newton–Raphson update for logistic regression, and it is also known as iteratively reweighted least squares (IRLS) because each step solves a weighted least squares problem with weights that change from iteration to iteration.

The practical difference from gradient ascent is striking. Newton–Raphson typically converges in four to eight iterations, regardless of the number of features, and requires no tuning of a learning rate. The price is the computation of (XᵀWX)⁻¹, which is an p × p matrix inversion performed once per iteration. For problems with thousands of features this becomes expensive, which is why quasi-Newton methods like L-BFGS — which approximate the Hessian without forming it — are the default in most libraries.

3.9 Why the maximum is unique: convexity

Both gradient ascent and Newton–Raphson find a maximum. Do they find the maximum — the global one — or could there be multiple peaks, with the algorithm getting stuck in a lesser one?

For logistic regression, the answer is guaranteed. The Hessian is −XᵀWX, and this matrix is always negative semi-definite, because:

• W is diagonal with positive entries p_i(1 − p_i) > 0.
• For any vector v, vᵀ XᵀWX v = (Xv)ᵀW(Xv) = Σ w_i (x_iᵀv)² ≥ 0.

A function whose Hessian is negative semi-definite everywhere is concave. A concave function has a single peak (or a flat plateau), with no local maxima and no saddle points. This means:

• Any local maximum is a global maximum.
• Gradient ascent will converge to the global maximum for a suitably small step size.
• Newton–Raphson, once it takes a step, converges rapidly and reliably.
• Different initialisations lead to the same answer (up to numerical precision).

This is a remarkable property, and it is not shared by all machine learning models. Neural networks, for example, have highly non-convex loss surfaces with many local minima. The convexity of logistic regression is what makes it such a reliable and well-understood model, and it is the reason the fitting procedure can be trusted to find the best answer.

3.10 Practical notes: fitting in scikit-learn

In practice you will rarely implement the optimiser yourself. The LogisticRegression class in scikit-learn provides a range of solvers, and choosing between them is mostly a matter of problem size and structure.

Two defaults are worth knowing. First, LogisticRegression applies L2 regularisation by default (with strength controlled by C = 1.0), so the coefficients are shrunk slightly relative to the pure maximum-likelihood estimate. To obtain the unregularised MLE, set penalty=None. Second, features are not scaled automatically, so if your features differ substantially in magnitude you should standardise them first — either with StandardScaler or by using a pipeline.

If you see a ConvergenceWarning, the usual causes are unscaled features, a small max_iter, or (less commonly) perfectly separable data. Scaling almost always fixes the first two; the third is addressed by regularisation.

04 Visual illustration

Two figures make the mathematics concrete. The first shows how the log loss penalises individual predictions; the second shows the shape of the log-likelihood surface for a small dataset.

Figure 1 explains why logistic regression punishes confident mistakes. A prediction of p = 0.9 when y = 0 incurs a loss of −log(0.1) ≈ 2.30, while a prediction of p = 0.99 incurs −log(0.01) ≈ 4.61 — twice as much, even though the two predictions differ by only 0.09 in probability. The loss grows without bound as p approaches the wrong extreme.

Figure 2 shows the log-likelihood surface for a small dataset, as a function of a single coefficient. The curve is concave — one smooth peak, no other local maxima — which is the property that guarantees the fitting procedure succeeds.

Two features of Figure 2 are worth noting. The curve is smooth — no kinks, no discontinuities — because both the sigmoid and the logarithm are infinitely differentiable. And it has exactly one peak, which is the convexity result of Section 3.9 made visible. A numerical optimiser started anywhere on this curve will climb to the same maximum.

05 Worked example

We work through a small dataset by hand to make the mechanics concrete. The dataset has four observations and a single feature:

The labels alternate with the feature, so the data is not separable — no threshold on x perfectly separates the classes. The model is z = β₀ + β₁x, with p = σ(z).

ℓ(0, 0) = Σ [ y_i log(0.5) + (1 − y_i) log(0.5) ]
        = 4 × log(0.5)
        = 4 × (−0.6931)
        = −2.7726

∂ℓ/∂β₀ = Σ (y_i − p_i)
       = (0 − 0.5) + (1 − 0.5) + (0 − 0.5) + (1 − 0.5)
       = −0.5 + 0.5 − 0.5 + 0.5
       = 0

∂ℓ/∂β₁ = Σ (y_i − p_i) x_i
       = (−0.5)(−2) + (0.5)(−1) + (−0.5)(1) + (0.5)(2)
       = 1 − 0.5 − 0.5 + 1
       = 1.0

β₀ ← 0 + 0.5 × 0  = 0
β₁ ← 0 + 0.5 × 1  = 0.5

∂ℓ/∂β₁ = (−0.2689)(−2) + (1 − 0.3775)(−1)
       + (−0.6225)(1) + (1 − 0.7311)(2)

     = 0.5378 − 0.6225 − 0.6225 + 0.5378

     = −0.1694

β₁ ← 0.5 + 0.5 × (−0.1694) = 0.4153

H₁₁ = − Σ w_i x_i²
    = − 0.25 × ((−2)² + (−1)² + 1² + 2²)
    = − 0.25 × 10
    = −2.5

β₁ ← 0 − (1.0) / (−2.5) = 0.4

06 Interactive activity

import numpy as np

def sigmoid(z):
    # numerically stable for large |z|
    return np.where(z >= 0,
                    1 / (1 + np.exp(-z)),
                    np.exp(z) / (1 + np.exp(z)))

def log_likelihood(X, y, beta):
    p = sigmoid(X @ beta)
    # clip to avoid log(0)
    p = np.clip(p, 1e-12, 1 - 1e-12)
    return np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))

def gradient(X, y, beta):
    p = sigmoid(X @ beta)
    return X.T @ (y - p)

# Gradient ascent
def fit_gradient_ascent(X, y, alpha=0.1, n_iter=5000, tol=1e-8):
    beta = np.zeros(X.shape[1])
    ll_history = [log_likelihood(X, y, beta)]
    for _ in range(n_iter):
        beta = beta + alpha * gradient(X, y, beta)
        ll = log_likelihood(X, y, beta)
        ll_history.append(ll)
        if abs(ll_history[-1] - ll_history[-2]) < tol:
            break
    return beta, ll_history

# Compare with scipy
from scipy.optimize import minimize

def neg_ll(beta, X, y):
    return -log_likelihood(X, y, beta)

def neg_grad(beta, X, y):
    return -gradient(X, y, beta)

result = minimize(neg_ll, np.zeros(X.shape[1]), args=(X, y),
                  jac=neg_grad, method='BFGS')
beta_scipy = result.x

Three things to notice as you run this. First, gradient ascent with a fixed learning rate is sensitive to the choice of α; the same value that converges quickly on one dataset may oscillate or crawl on another. Second, scipy's BFGS reaches the same answer in far fewer iterations, because it approximates the curvature of the loss surface. Third, the coefficients from your implementation and from scikit-learn will match closely only if you remember to disable regularisation — by default, scikit-learn shrinks the coefficients.

07 Self-check questions

Answer each question before revealing the solution.

08 Practice exercises

Work these by hand where possible, then verify computationally.

1. For the worked example dataset, compute the log-likelihood at β₀ = 0, β₁ = 0.4 and at β₀ = 0, β₁ = 0.5. Which is higher? Does this match Figure 2?
2. Compute the gradient ∂ℓ/∂β₁ at β₁ = 0.4 for the worked example. Is it positive or negative? What does the sign tell you about whether to increase or decrease β₁?
3. A single observation has y = 1 and p = 0.2. What is the log loss? What would it be if p = 0.9?
4. Compute the Hessian entry H₁₁ for the worked example at β₁ = 0.4. Is it negative? What does this confirm?
5. Starting from β₁ = 0.4, take one Newton–Raphson step using the Hessian from exercise 4. How close does this get to the optimum of 0.42?
6. A logistic regression is fitted to a dataset with 200 observations and 3 features. What is the dimension of the Hessian matrix? How many entries does it have?
7. Show algebraically that σ′(z) = σ(z)(1 − σ(z)). This is the fact used in deriving the Hessian.
8. Consider a dataset where every label is 1. What is the maximum-likelihood estimate of the intercept for an intercept-only model? What happens to the coefficients on the features?

09 Problem sets

These problems require reasoning about the fitting procedure rather than calculation.

10 Discussion and projects

Discussion questions

1. Maximum likelihood is a general principle, not specific to logistic regression. What other models have you met that can be fitted by maximum likelihood? What is the same and what is different?
2. The convexity of the log-likelihood means logistic regression has a unique best answer. Does this make it a "better" model than one with a non-convex loss surface, such as a neural network? Or does it make it a more limited one?
3. The log loss penalises confident mistakes without bound. Is this always the right behaviour? Can you think of a context in which penalising confidence too harshly would be undesirable?
4. In practice, almost nobody uses unregularised logistic regression. Given that the maximum-likelihood estimate is theoretically optimal (in the sense of asymptotic efficiency), why is the regularised version preferred?

Mini-project · Implement and compare optimisers

In this project you will implement three different ways of fitting logistic regression and compare them empirically. The goal is not to build a production tool but to understand what the library does on your behalf.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler

# Generate and scale data
X, y = make_classification(n_samples=400, n_features=4, n_informative=3,
                           n_redundant=0, random_state=0)
X = StandardScaler().fit_transform(X)
X = np.hstack([np.ones((X.shape[0], 1)), X])  # add intercept column

def sigmoid(z):
    return np.where(z >= 0,
                    1 / (1 + np.exp(-z)),
                    np.exp(z) / (1 + np.exp(z)))

def log_likelihood(beta, X, y):
    p = np.clip(sigmoid(X @ beta), 1e-12, 1 - 1e-12)
    return np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))

def gradient(beta, X, y):
    return X.T @ (y - sigmoid(X @ beta))

def hessian(beta, X, y):
    p = sigmoid(X @ beta)
    W = np.diag(p * (1 - p))
    return -X.T @ W @ X

# Gradient ascent
def gradient_ascent(X, y, alpha=0.1, n_iter=500):
    beta = np.zeros(X.shape[1])
    history = [log_likelihood(beta, X, y)]
    for _ in range(n_iter):
        beta = beta + alpha * gradient(beta, X, y)
        history.append(log_likelihood(beta, X, y))
    return beta, history

# Newton-Raphson
def newton_raphson(X, y, n_iter=20):
    beta = np.zeros(X.shape[1])
    history = [log_likelihood(beta, X, y)]
    for _ in range(n_iter):
        g = gradient(beta, X, y)
        H = hessian(beta, X, y)
        beta = beta - np.linalg.solve(H, g)
        history.append(log_likelihood(beta, X, y))
    return beta, history

beta_ga, history_ga = gradient_ascent(X, y)
beta_nr, history_nr = newton_raphson(X, y)

fig, ax = plt.subplots(figsize=(8, 4.8))
ax.plot(history_ga, label="gradient ascent", linewidth=2)
ax.plot(history_nr, label="Newton-Raphson", linewidth=2)
ax.set_xlabel("iteration")
ax.set_ylabel("log-likelihood")
ax.set_title("Convergence comparison")
ax.legend()
ax.grid(alpha=0.25)
plt.tight_layout()
plt.show()

11 Summary

This tutorial answered the question of where the coefficients come from.

• Least squares is the wrong criterion. It under-penalises confident mistakes and does not correspond to the Bernoulli probability model assumed by logistic regression.
• Maximum likelihood is the right criterion. Choose the coefficients that make the observed labels most probable under the model.
• The likelihood for a binary observation is p^y (1 − p)^(1−y). Taking logs turns the product over observations into a sum, giving the log-likelihood.
• The negative log-likelihood is the log loss. It penalises confident mistakes without bound, which is exactly the behaviour a loss function for probabilities should have.
• The gradient is Xᵀ(y − p). Each component is a weighted sum of residuals — the same structural form as in linear regression, with fitted probabilities replacing fitted values.
• There is no closed-form solution. The normal equations are non-linear in β because of the sigmoid, so the coefficients must be found iteratively.
• Gradient ascent works but is slow. It requires a suitable learning rate, and its speed depends on the scaling of the features.
• Newton–Raphson is much faster. Using the exact Hessian −XᵀWX, it converges in a handful of iterations and needs no learning rate.
• The log-likelihood is concave. The Hessian is negative semi-definite everywhere, so any local maximum is the global maximum, and the fitted coefficients are unique.
• scikit-learn handles all of this. It provides several solvers, applies L2 regularisation by default, and does not scale your features for you.

12 Further reading

Core

• James, Witten, Hastie & Tibshirani — An Introduction to Statistical Learning, Section 4.3. Covers logistic regression and maximum likelihood estimation, with derivations of the log-likelihood and its properties.
• Hastie, Tibshirani & Friedman — The Elements of Statistical Learning, Section 4.4. A more rigorous treatment including the Newton–Raphson algorithm, IRLS, and the connection to generalised linear models.
• scikit-learn documentation — LogisticRegression. Read the "solvers" section for the practical trade-offs between L-BFGS, Newton methods, coordinate descent and stochastic gradient methods.

Supplementary

• Bishop — Pattern Recognition and Machine Learning, Section 4.3. Develops logistic regression within a probabilistic framework, with a careful derivation of the log-likelihood and the IRLS algorithm.
• Murphy — Probabilistic Machine Learning: An Introduction, Chapter 8. A modern treatment with clear discussion of both the mathematical foundations and the practical optimisation issues.
• McCullagh & Nelder — Generalized Linear Models. The classic reference. Logistic regression is a special case of a GLM with a binomial response and a logit link, and this book develops the general theory that unifies it with linear regression and other models.

On optimisation

• Boyd & Vandenberghe — Convex Optimization. The standard reference for the convexity arguments used in Section 3.9. Chapters 2–4 establish the properties of convex and concave functions and their implications for optimisation.
• Nocedal & Wright — Numerical Optimization. A comprehensive treatment of gradient methods, Newton methods, and quasi-Newton methods such as BFGS and L-BFGS.
• King & Zeng — "Logistic Regression in Rare Events Data" (2001). A well-known paper on the practical difficulties of fitting logistic regression when the positive class is rare — a preview of the material in T-06.