Generalized Linear Models (GLM)
The exponential family, link functions, and IRLS as Newton's method in disguise.
On this page
Beginner: ordinary least squares (2.3.1) quietly assumes something specific about the target y: that it's Gaussian-distributed, scattered symmetrically around a mean, and that mean happens to be a straight-line (linear) function of the features. That's a fine assumption for, say, predicting a house's price in thousands of dollars. It's a terrible assumption the moment y is a count — the number of insurance claims filed this year, which can never go negative and is lumpy, not smoothly bell-shaped — or a probability of a binary outcome, like whether a patient has a disease, which must stay between 0 and 1 no matter what the features say. A generalized linear model (GLM) keeps the "linear combination of features" part everyone likes about OLS, but frees up the distributional assumption: y can come from any distribution in a large, well-behaved family called the exponential family (which happens to include Gaussian, Bernoulli, Poisson, Gamma, and several others), and a chosen link function connects the ordinary linear predictor Xθ to that distribution's mean.
Intermediate: the link function is the whole trick. Instead of directly setting the mean equal to Xθ (which breaks the moment the mean is constrained — a probability can't be 1.4, a count can't be −3), a GLM sets some transformed version of the mean equal to Xθ, and lets that transformation absorb the constraint. The identity link (no transformation at all) recovers plain OLS exactly. The logit link, log(p/(1-p)), maps a probability in (0,1) onto the whole real line, so a linear predictor of any sign is always mapped back to a valid probability — this is exactly logistic regression, covered in full in a later Classification module. The log link maps a positive count's mean onto the real line the same way, giving Poisson regression. Three seemingly unrelated models — linear, logistic, Poisson regression — are one model with three different link/distribution choices plugged in.
Advanced: what makes this more than a naming convenience is that every member of the exponential family shares the same algebraic skeleton (the Formula section below), and fitting any GLM by maximum likelihood reduces to the same algorithm: Iteratively Reweighted Least Squares (IRLS). The Derivation section below shows this isn't a coincidence or a separate technique that happens to resemble Newton's method — IRLS is Newton's method (derived in section 2.2.4), applied to the GLM log-likelihood and simplified until its Hessian is recognized as reducible to a weighted least-squares normal-equations solve. One optimization idea, one family of distributions, one link mechanism — and OLS, logistic regression, and Poisson regression fall out as three special cases rather than three separate models to memorize.
Worked intuition — three link functions, one framework. Suppose the same linear predictor η = Xθ comes out to 2.0 for some input. Under the identity link (Gaussian family), the predicted mean is literally 2.0 — a continuous, unconstrained quantity, so no transformation is needed. Under the logit link (Bernoulli family), that same 2.0 is passed through the logistic sigmoid to give a predicted probability of about 0.88 — always between 0 and 1, however extreme η gets. Under the log link (Poisson family), it's exponentiated to give a predicted mean count of about 7.4 — always positive. Same linear predictor, three different guarantees about what comes out the other side, purely from the choice of g.
The GLM link equation — a link function g connects the linear predictor to the mean of y:
The exponential-family density that every GLM assumes for y, in its canonical form (η is the natural parameter, φ a dispersion parameter, and a, b, c are known functions that vary by distribution):
Every distribution used below is this same shape with different a, b, c plugged in: Gaussian with the identity link, Bernoulli with the logit link, Poisson with the log link, Gamma with the inverse (or log) link. When the link is chosen so that η — the natural parameter above — equals the linear predictor Xθ directly, it's called the canonical link for that family, and it is exactly what makes the derivation below come out so cleanly (more on non-canonical links in the Expert note).
Setup. Take n independent observations under a canonical link, so the natural parameter for observation i is exactly the linear predictor, η_i = x_i^Tθ. The dispersion term a(φ) is a positive scalar shared by every observation and doesn't depend on θ at all, and c(y,φ) doesn't depend on θ either — so for the purposes of differentiating with respect to θ, both can be dropped (a(φ) will be seen to cancel exactly in Newton's update below anyway). What's left of the log-likelihood, summed over the sample:
(a) The gradient. Differentiate term i with respect to η_i first, then apply the chain rule through η_i = x_i^Tθ (so ∂η_i/∂θ = x_i):
Define μ_i = b'(η_i) — this is exactly the mean of y_i under the exponential-family density above, a standard identity for this family (the same way a Gaussian's mean is its own natural parameter). Written as one vector equation over the whole design matrix X:
— precisely the "prediction error times features" gradient form that shows up for OLS and logistic regression separately; here it's shown to be the same formula for the entire exponential family at once.
(b) The Hessian. Differentiate the gradient once more with respect to θ. Each term's dependence on θ runs entirely through η_i again:
where W is diagonal with entries w_i = b''(η_i) — the variance-related weight for each observation (for the exponential family, b''(η_i) is exactly the variance of y_i, up to the shared dispersion scalar dropped above). Every entry of W depends on the current θ through η_i — which is exactly why this has to be an iterative, re-weighted procedure rather than a one-shot linear solve.
(c) Substitute into Newton's update. Newton's method's update rule, derived from a second-order Taylor expansion in section 2.2.4, is:
Plugging in the gradient and Hessian just derived:
This is already a Newton step, but it doesn't yet look like a least-squares problem. To see that it is one, left-multiply both sides by X^TWX and regroup:
(the last step just factors X^TW out, using X^T(y-μ) = X^TW·W^{-1}(y-μ), which is trivially true since W·W^{-1} = I). Naming that bracketed quantity the working response:
the Newton update collapses to:
Compare this to section 2.3.1's ordinary normal equations, θ = (X^TX)^{-1}X^Ty. The IRLS update is exactly that formula with two substitutions: an observation-weighted design (W inserted between every X^T and X) and a working response z in place of the raw target y. That is the entire algorithm: recompute μ, W, and z from the current θ, solve a weighted OLS problem for a new θ, and repeat until θ stops moving. "Iteratively Reweighted Least Squares" is a literal description of what the loop does, not a separate algorithm that happens to resemble Newton's method — it is Newton's method, with its particular Hessian recognized as a weighted normal-equations solve.
The Gaussian/identity-link case converges in exactly one step. For the Gaussian family under its canonical (identity) link, b(η) = ½η², so μ = b'(η) = η and w = b''(η) = 1 for every observation — a constant weight, so W = I regardless of θ. Because the link is the identity, μ = η = Xθ_{old} exactly, so the working response simplifies immediately:
The working response is just y itself, no matter what θ_old was. Substituting W = I and z = y into the IRLS update:
— the plain OLS normal equations, arrived at on the very first iteration, from any starting θ_old (even θ_old = 0), with zero iteration needed beyond that. The "generalized" machinery hasn't just approximated OLS — it has degenerated back into being OLS exactly, because a Gaussian/identity-link GLM's log-likelihood is already a plain quadratic in θ, and Newton's method solves a quadratic exactly in one step (the same fact section 2.2.4 uses to explain why Newton's method converges so fast near any well-behaved minimum). Every non-Gaussian or non-identity-link case is where W and z genuinely change from round to round, and the "iteratively" in IRLS starts to earn its name.
Where this is used: when you call R's glm(), Python's statsmodels.api.GLM(...).fit(), or almost any other statistical package's GLM routine, the loop it runs internally — by default — is precisely this one: recompute weights and a working response from the current fit, solve a weighted least-squares problem, repeat. Everything above is not a simplified teaching version of what production statistical software does; it is what production statistical software does.
Six IRLS/Newton rounds, replayed automatically, fitting a logistic-link GLM to 14 overlapping binary-labeled points. Each round re-solves the weighted least-squares step derived above and redraws the fitted probability curve; the log-likelihood readout below climbs every round on this well-behaved, non-separable dataset -- see the note under the diagram for why that is not the same guarantee EM's ELBO argument makes.
The Worked-intuition example above, made interactive: drag η and pick a link to watch the same linear predictor pass through g⁻¹ into a very differently-shaped response — an unbounded line for the identity link, a bounded S-curve for the logit link, and an always-positive exponential for the log link.
All three tabs fit the same canonical-link Bernoulli GLM (logistic regression) to the same 14-point dataset used in the diagram above. The first two hand-roll IRLS exactly as derived — recompute weights and a working response, solve the weighted normal equations, repeat — with no logistic-regression or GLM library anywhere. The third swaps in statsmodels' real GLM class (which runs the identical IRLS loop internally) to confirm the from-scratch coefficients match a production implementation, then fits a second, log-link Poisson GLM on synthetic count data to show the same machinery handling a different exponential-family member.
#include <cmath>
#include <cstdio>
#include <vector>
#include <array>
// Same IRLS loop as the Python version, with the 2x2 weighted normal equations
// solved directly by hand (no linear-algebra library at all).
static const std::vector<double> xs = {-4, -3, -2, -1.5, -1, -0.5, 0,
0.5, 1, 1.5, 2, 2.5, 3, 4};
static const std::vector<double> ys = {0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1};
double sigmoid(double eta) { return 1.0 / (1.0 + std::exp(-eta)); }
double logLikelihood(const std::array<double, 2>& theta) {
double s = 0.0;
for (size_t i = 0; i < xs.size(); ++i) {
double eta = theta[0] + theta[1] * xs[i];
double log1pExpEta = eta > 30.0 ? eta : std::log1p(std::exp(eta));
s += ys[i] * eta - log1pExpEta;
}
return s;
}
// One IRLS round: build the 2x2 weighted normal equations (X^T W X) theta = X^T W z
// and solve them with Cramer's rule -- an explicit weighted-least-squares solve,
// exactly the update derived in the lesson.
std::array<double, 2> irlsStep(const std::array<double, 2>& theta) {
double a00 = 0, a01 = 0, a11 = 0, b0 = 0, b1 = 0;
for (size_t i = 0; i < xs.size(); ++i) {
double x = xs[i];
double eta = theta[0] + theta[1] * x;
double mu = sigmoid(eta);
double w = std::max(mu * (1.0 - mu), 1e-6);
double z = eta + (ys[i] - mu) / w;
a00 += w;
a01 += w * x;
a11 += w * x * x;
b0 += w * z;
b1 += w * x * z;
}
double det = a00 * a11 - a01 * a01;
return {(a11 * b0 - a01 * b1) / det, (a00 * b1 - a01 * b0) / det};
}
int main() {
std::array<double, 2> theta = {0.0, 0.0};
std::printf("round 0: theta=(%.4f, %.4f) log-likelihood=%.4f\n",
theta[0], theta[1], logLikelihood(theta));
for (int round = 1; round <= 7; ++round) {
theta = irlsStep(theta);
std::printf("round %d: theta=(%.4f, %.4f) log-likelihood=%.4f\n",
round, theta[0], theta[1], logLikelihood(theta));
}
return 0;
}
// Converges to the same theta ~= (-0.348, 1.427) as the Python version, in the
// same handful of rounds -- the algorithm is entirely deterministic linear algebra,
// so language choice changes nothing about the trajectory.- Logistic regression (logit link, Bernoulli family) for binary outcomes — medical diagnosis (disease present or not, given clinical measurements) and churn/attrition prediction (will this customer cancel this month) are the flagship real-world use of a GLM. Its full treatment as a standalone classifier — decision boundaries, calibration, evaluation metrics — is developed in a later Classification module; here it's the concrete example that motivated the entire IRLS derivation above.
- Poisson regression (log link, Poisson family) for count data — the number of insurance claims a policyholder files in a year, the number of visits/events a web page receives per day, or the number of new cases of a disease reported per region per week. All three are non-negative integer counts with no natural upper bound, exactly the shape the Poisson family is built for and plain Gaussian OLS is not.
- Gamma-family GLMs (typically a log or inverse link) for claim severity — insurance pricing usually splits into two separate GLMs: a Poisson (or negative-binomial) GLM for how many claims occur (frequency), and a Gamma GLM for how expensive each claim is, given that it occurred (severity). Claim amounts are strictly positive and typically right-skewed (many small claims, a long tail of large ones) — exactly what the Gamma family models, and exactly what a symmetric Gaussian assumption gets wrong.
- The identity-link Gaussian GLM is exactly ordinary least squares. Section 2.3.1's OLS regression isn't a different technique sitting next to GLMs — it's the special case worked out explicitly in the Derivation above, where
Wis constant and IRLS converges in exactly one step. Every intuition already built about OLS carries over unchanged; it's just now visible as one corner of a much larger framework. - Multinomial/softmax regression extends the same idea to more than two classes (a categorical, rather than Bernoulli, response) — one linear predictor per class, tied together through a softmax link. It's covered separately as a classification model in a later module, but structurally it's one more exponential- family member with a matching canonical link, following the same pattern.
- Every standard statistics package's
glm()/GLM()call — R's built-inglm(), Python'sstatsmodels.api.GLM, SAS'sPROC GENMOD, and others all expose one function that takes a family and a link as arguments and runs the IRLS loop derived above underneath, regardless of which of the above examples it's being asked to fit.
- Fitting plain OLS to count, probability, or strictly-positive data anyway. This is a very common real-world mistake, and it produces genuinely broken predictions, not just "slightly suboptimal" ones: an OLS model fit to a binary 0/1 outcome will happily predict probabilities below 0 or above 1 for some inputs, and an OLS model fit to count data can predict a negative expected claim count — both nonsensical for the quantity being modeled, and both symptoms of forcing a Gaussian/identity-link assumption onto data that was never Gaussian to begin with. The fix is exactly this lesson: pick the family that actually matches the response's support and shape, and let the link function handle keeping predictions in range.
- Confusing the link function with the distributional family. They are two independent choices, not one — you could pair a Bernoulli family with a probit link instead of the logit link, for instance, and get a different (but perfectly valid) model. Canonical links, as derived above, are simply the choice that makes
ηcoincide with the family's own natural parameter, which is what gives IRLS its cleanest algebra (see the Expert note below for what changes with a non-canonical link) — "canonical" describes a convenient algebraic property, not the only valid or even always the best choice for a given family. - Assuming IRLS always converges in a single step, the way the Gaussian/identity case does. That one-step convergence is a special property of a constant
W, not a general feature of IRLS — for logistic, Poisson, Gamma, and every other non-Gaussian family,Wand the working response genuinely change every round, and several iterations are typically needed. Worse, on linearly separable data — a real, well-known failure mode of logistic regression specifically — the fitted probabilities can be pushed arbitrarily close to 0 or 1 without ever finishing, the weightsw_i = μ_i(1-μ_i)collapse toward 0, the weighted normal equations become numerically singular, and the coefficients diverge toward infinity instead of converging at all.
Going deeper
Canonical vs. non-canonical links, precisely. A link g is canonical for a given exponential-family member exactly when it makes the natural parameter equal to the linear predictor, η = Xθ, with no further transformation in between — which is precisely the assumption the Derivation section above leaned on to get such a clean gradient (X^T(y-μ)) and Hessian (-X^TWX). Choose a non-canonical link instead — say, a probit link for a Bernoulli outcome instead of the logit — and μ = g^{-1}(η) no longer coincides with b'(η) for the family's own natural parameter. The chain rule then has to route through an extra factor, dμ/dη, at both the gradient and the Hessian stage, and the exact (observed) Hessian picks up an additional second-derivative term that the canonical-link derivation above never had to deal with — the weights are no longer simply b''(η).
This is precisely why most production GLM implementations don't literally run raw Newton-Raphson on the observed Hessian for non-canonical links — they run Fisher scoring instead, which replaces the observed Hessian with its expected value under the model. For a canonical link, the observed and expected information coincide exactly (which is exactly why the derivation above needed no such distinction), so Newton-Raphson and Fisher scoring are the same algorithm. For a non-canonical link, they differ, and Fisher scoring is generally preferred because the resulting weights are guaranteed non-negative and the iteration tends to be numerically steadier — a genuine wrinkle worth knowing before reading a real GLM library's source.
For an identity-link Gaussian GLM, IRLS converges in exactly one iteration from any starting theta. Why does that specific shortcut not carry over to a logit-link Bernoulli GLM (logistic regression), even though both are fit by the same IRLS algorithm?
The one-step shortcut relies on two facts that are both true only for the Gaussian/identity case: the weight w_i = b''(eta_i) = 1 is a constant, so W = I regardless of theta, and because the link is the identity, mu = eta = X*theta_old exactly, so the working response z = X*theta_old + (y - X*theta_old) collapses to just y itself -- independent of theta_old. That means round one's weighted least-squares problem is already the plain, unweighted OLS problem, so its solution is the final answer immediately. For a logit-link Bernoulli GLM, w_i = mu_i*(1-mu_i) depends on the current fit's predicted probability and is not constant, and mu is a nonlinear (sigmoid) function of eta rather than equal to it, so z = eta_old + (y-mu)/w does not simplify to a fixed target independent of theta_old. Each round genuinely produces a different weighted least-squares problem with different weights and a different working response, so several rounds are generally needed before theta stops moving -- exactly the 'iteratively' in Iteratively Reweighted Least Squares earning its name.
A GLM keeps regression's "linear combination of features" backbone but frees up two choices that OLS bakes in as fixed: which exponential-family distribution generates y, and which link function connects the linear predictor to that distribution's mean. Fitting any GLM by maximum likelihood reduces to one algorithm, IRLS, which the Derivation section above showed is exactly Newton's method (2.2.4) applied to the GLM log-likelihood — its Hessian recognized as reducible to a weighted least-squares normal-equations solve, repeated against a working response that updates every round. The identity-link Gaussian case collapses that loop to a single step and recovers OLS (2.3.1) exactly; the logit-link Bernoulli case is logistic regression; the log-link Poisson case is Poisson regression — three named models that turn out to be one model wearing three different link/family choices.