Generalized linear models (GLMs)

Deep dive

How it works

Parts:

• Family — distribution of y, mean μ, variance φ·V(μ): Poisson μ, gamma μ², binomial μ(1−μ), Tweedie μ^q (1 < q < 2).
• Link g — g(μ) = η; log for counts and amounts, logit for proportions. Canonical link: Xᵀy is sufficient for β.
• Linear predictor and offset — η = Xβ + offset; the offset is a known per-row log(exposure) (policy-years, person-time, population) with coefficient 1.
• IRLS loop — working response and weights, then the weighted normal equations by Cholesky.
• Deviance — 2·(saturated log-likelihood − fitted log-likelihood); goodness of fit and stopping criterion.

Operations (one IRLS iteration, current η):

1. μ = g⁻¹(η).
2. Working weight w = 1 / (V(μ)·g'(μ)²); working response z = η − offset + (y − μ)·g'(μ).
3. Solve (XᵀWX) β = XᵀWz.
4. η = Xβ + offset; compute deviance; stop when its relative change is below about 1e-8.

Start from μ0 = y + small so η0 = g(μ0) is finite. Afterwards Cov(β̂) = φ̂·(XᵀWX)⁻¹, with φ = 1 for Poisson and binomial, else φ̂ = Pearson χ² / (n − p). Weights: Poisson with log link w = μ; gamma with log link w = 1; binomial with logit w = μ(1−μ).

 y, X, offset ──► η = Xβ+off ──► μ = g⁻¹(η) ──► w, z ──► solve XᵀWX β = XᵀWz ─┐
                      ▲                                                          │
                      └────────────────── new β ◄── deviance change < tol? ◄────┘

Code

Compiled and run; all asserts pass. Synthetic data, known β, exposure correlated with x1.

/// Poisson, log link, offset, IRLS. x: n×p row-major, column 0 = intercept.
fn poisson_irls(x: &[f64], y: &[f64], off: &[f64], p: usize) -> (Vec<f64>, Vec<f64>) {
    let n = y.len();
    let mut eta: Vec<f64> = y.iter().map(|&v| (v + 0.1).ln()).collect(); // mu0 = y + 0.1
    let (mut beta, mut devs) = (vec![0.0; p], Vec::<f64>::new());
    for _ in 0..25 {
        let (mut a, mut b) = (vec![0.0; p * p], vec![0.0; p]);
        for i in 0..n {
            let mu = eta[i].exp();
            let z = eta[i] - off[i] + (y[i] - mu) / mu; // working response (offset removed)
            let r = &x[i * p..(i + 1) * p];
            for j in 0..p { // weight w = 1/(V g'^2) = mu
                b[j] += mu * r[j] * z;
                for k in 0..=j { a[j * p + k] += mu * r[j] * r[k]; }
            }
        }
        beta = chol_solve(a, &b, p);
        let mut dev = 0.0;
        for i in 0..n {
            let r = &x[i * p..(i + 1) * p];
            eta[i] = off[i] + (0..p).map(|j| r[j] * beta[j]).sum::<f64>();
            let mu = eta[i].exp();
            dev += 2.0 * (if y[i] > 0.0 { y[i] * (y[i] / mu).ln() } else { 0.0 } - (y[i] - mu));
        }
        let done = devs.last().map_or(false, |d| (d - dev).abs() < 1e-8 * (dev + 0.1));
        devs.push(dev);
        if done { break; }
    }
    (beta, devs)
}

/// Solve A z = b by Cholesky (A symmetric positive definite, lower triangle filled).
fn chol_solve(mut a: Vec<f64>, b: &[f64], p: usize) -> Vec<f64> {
    for j in 0..p {
        let s = a[j * p + j] - (0..j).map(|k| a[j * p + k].powi(2)).sum::<f64>();
        assert!(s > 0.0, "X'WX not positive definite");
        a[j * p + j] = s.sqrt();
        for i in j + 1..p {
            let t = a[i * p + j] - (0..j).map(|k| a[i * p + k] * a[j * p + k]).sum::<f64>();
            a[i * p + j] = t / a[j * p + j];
        }
    }
    let mut z = b.to_vec();
    for i in 0..p { for k in 0..i { z[i] -= a[i * p + k] * z[k]; } z[i] /= a[i * p + i]; }
    for i in (0..p).rev() { for k in i + 1..p { z[i] -= a[k * p + i] * z[k]; } z[i] /= a[i * p + i]; }
    z
}

fn main() {
    let (n, p, truth) = (20_000, 4, [-0.5, 0.5, -0.3, 0.2]); // truth: log rate per unit exposure
    let mut s = 0x2545F4914F6CDD1Du64;
    let mut u = move || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s >> 11) as f64 / (1u64 << 53) as f64 };
    let (mut x, mut y, mut off) = (vec![1.0; n * p], vec![0.0; n], vec![0.0; n]);
    for i in 0..n {
        for j in 1..p { x[i * p + j] = 2.0 * u() - 1.0; }
        off[i] = 0.7 * x[i * p + 1] + 0.3 * u(); // log exposure, correlated with x1
        let eta: f64 = off[i] + (0..p).map(|j| x[i * p + j] * truth[j]).sum::<f64>();
        let (l, mut t) = ((-eta.exp()).exp(), 1.0); // Knuth Poisson sampler
        while { t *= u(); t > l } { y[i] += 1.0; }
    }
    let (b, devs) = poisson_irls(&x, &y, &off, p);
    println!("beta = {b:.3?}\ndeviance per pass = {devs:.2?}");
    assert!(b.iter().zip(&truth).all(|(a, t)| (a - t).abs() < 0.05));
    assert!(devs.windows(2).all(|w| w[1] <= w[0] + 1e-6 * w[0]) && devs.len() <= 8);
    // doubling every exposure: slopes unchanged, intercept moves by -ln 2
    let off2: Vec<f64> = off.iter().map(|o| o + 2f64.ln()).collect();
    let (b2, _) = poisson_irls(&x, &y, &off2, p);
    assert!((b2[0] - (b[0] - 2f64.ln())).abs() < 1e-6 && (1..p).all(|j| (b2[j] - b[j]).abs() < 1e-6));
    // dropping the offset leaks the exposure effect into beta1 (0.5 + ~0.7)
    let (b0, _) = poisson_irls(&x, &y, &vec![0.0; n], p);
    println!("no offset: beta = {b0:.3?}");
    assert!(b0[1] - b[1] > 0.5);
}

Output: beta = [-0.497, 0.487, -0.303, 0.197] (truth [-0.5, 0.5, -0.3, 0.2]); deviance per pass [24702.30, 20752.89, 20556.05, 20554.91, 20554.91] (5 passes, the last confirms); without offset β1 = 1.186: the exposure effect (0.7 on x1) leaks in.

Complexity and constants

• Per iteration: n·p(p+1)/2 multiply-adds for XᵀWX, one exp per row, p³/3 for Cholesky; for n ≫ p the pass over rows dominates.
• Measured on this Mac (Apple M4 Pro), one core, release build, code above at n = 1M, p = 10: 27 ms per IRLS pass (about 27 ns per row), 5 passes, 0.135 s total, max coefficient error 0.003. The loop is scalar and untuned.

Efficiency

• Cost: about 27 ns per row per IRLS pass at p = 10 (measured on this Mac, M4 Pro, one core): one exp, 55 multiply-adds for the lower triangle of XᵀWX, a deviance pass. The p³/3 solve is negligible for p ≪ n. Bound: the p² accumulation (compute).
• Levers:
  • Aggregate identical rows (rating cells) — cost falls from n to the number of cells; exact for Poisson with offset, as y and exposure add (math) — categorical factors.
  • Newton-Cholesky instead of L-BFGS — scikit-learn docs call it a good choice for n_samples ≫ n_features; lbfgs stays default "for its robustness" — dense, modest p.
  • Parallel row chunks — independent partial XᵀWX sums, scaling with cores until memory bandwidth binds (not measured here) — large n.
  • Warm-started λ path by coordinate descent — whole path for about the cost of one fit (Friedman et al. 2010) — L1/elastic net.
  • FTRL or L-BFGS on sparse features — O(nnz) per pass, no p×p matrix — very large p.
• Naive → standard → state of the art:
  • Explicit (XᵀWX)⁻¹ each iteration · O(n·p² + p³), unstable with collinear columns
  • IRLS with Cholesky or QR, aggregated rows, step halving · few passes, O(n·p²) each
  • Multithreaded chunked accumulation, path-wise sparse solvers (glmnet, glum), L-BFGS for huge p · O(nnz) per pass

Implementation notes

• Cost drivers: one pass over the data per iteration, the p² accumulation, one exp per row. Accumulate in f64 (f32 sums over millions of rows can lose positive definiteness), in row chunks that fit L2 cache.
• Clamp μ away from 0 (and 1 for binomial) and η to about ±30 (log link) before computing w and z.
• Starting values: μ0 = y + 0.1 (Poisson); (y + 0.5)/2 for binomial 0/1.
• Non-canonical links and bad starts can raise the deviance: halve the step until it falls.
• Offset (count y, offset log e) and weight (rate y/e, weight e) give the same Poisson estimates.

Pitfalls

• Forgetting the offset biases every covariate correlated with exposure (run above: β1 0.49 → 1.19). Doubling all exposures leaves slopes unchanged and moves the intercept by −ln 2 (asserted).
• Overdispersion: variance above the mean makes Poisson standard errors too small. Use quasi-Poisson (scale by φ̂), negative binomial (V = μ + μ²/θ, a GLM for fixed θ) or Tweedie.
• Separation in logistic regression: the MLE does not exist and coefficients diverge. Add an L2 penalty (regularization).
• Training deviance is not accuracy: check out-of-sample deviance and predicted over observed per bin (calibration).
• Tweedie power is not fitted by IRLS: pick it by profile likelihood or cross-validation.

Variants and alternatives

In production

Actuarial pricing, incidence models, ecology counts, ad click and conversion models, defect-rate monitoring and credit scoring, fitted with R glm, statsmodels, scikit-learn, glum, glmnet and H2O.

Implementations

Sources

• Nelder and Wedderburn, Generalized Linear Models, Journal of the Royal Statistical Society A 135(3), 1972. paper
• McCullagh and Nelder, Generalized Linear Models, 2nd ed., Chapman & Hall 1989.
• Friedman, Hastie, Tibshirani, Regularization Paths for Generalized Linear Models via Coordinate Descent, Journal of Statistical Software 2010. paper
• Jørgensen, Exponential Dispersion Models, Journal of the Royal Statistical Society B 49(2), 1987 (Tweedie family).