Linear and logistic regression

Deep dive

How it works

In: feature matrix X (n × d), targets y, penalty λ · Out: weights w and intercept b · State: none after fitting (w, b are the model) · Cost: O(n·d² + d³) closed form or per Newton step

Ridge (least squares + L2).

1. Center or add an intercept column of ones. Do not penalize the intercept.
2. Form A = XᵀX + λI (d×d, symmetric positive definite for λ > 0) and c = Xᵀy.
3. Cholesky A = LLᵀ, then forward and back substitution. d³/3 flops for the factorization.

Logistic regression by Newton / IRLS.

1. p = σ(Xw). Gradient g = Xᵀ(p − y) + λw. Hessian H = XᵀWX + λI, W_ii = p_i(1 − p_i).
2. Solve H·Δ = g by Cholesky, update w ← w − Δ.
3. Stop when ‖Δ‖ or the relative loss change is below about 1e-8. Quadratic convergence near the optimum: typically 5-20 iterations. Each step is a weighted least-squares solve, hence "iteratively reweighted least squares".

Large n or d. Mini-batch SGD or L-BFGS on the same gradients, or coordinate descent (liblinear, glmnet). Newton is O(d³) per step, fine to d ≈ order of 10⁴.

          linear                         logistic
 x ─► w·x+b ─► ŷ             x ─► w·x+b = z ─► σ(z) = p ─► log loss
 loss ½(ŷ−y)², grad (ŷ−y)x    loss −y ln p − (1−y) ln(1−p), grad (p−y)x

Code

/// Dense row-major helpers. X: n×d, A: d×d.
fn xtwx(x: &[f64], wts: &[f64], n: usize, d: usize, lam: f64) -> Vec<f64> {
    let mut a = vec![0.0; d * d];
    for i in 0..n {
        let r = &x[i * d..(i + 1) * d];
        for j in 0..d { let s = wts[i] * r[j]; for k in 0..=j { a[j * d + k] += s * r[k]; } }
    }
    for j in 0..d { for k in 0..j { a[k * d + j] = a[j * d + k]; } }
    for j in 1..d { a[j * d + j] += lam; }                     // column 0 = intercept: not penalized
    a
}

/// Solve A z = b for symmetric positive definite A via Cholesky (A = L Lᵀ).
fn cholesky_solve(mut a: Vec<f64>, b: &[f64], d: usize) -> Vec<f64> {
    for j in 0..d {
        let mut s = a[j * d + j];
        for k in 0..j { s -= a[j * d + k] * a[j * d + k]; }
        assert!(s > 0.0, "matrix not positive definite: add lambda");
        a[j * d + j] = s.sqrt();
        for i in j + 1..d {
            let mut t = a[i * d + j];
            for k in 0..j { t -= a[i * d + k] * a[j * d + k]; }
            a[i * d + j] = t / a[j * d + j];
        }
    }
    let mut z = b.to_vec();
    for i in 0..d { for k in 0..i { z[i] -= a[i * d + k] * z[k]; } z[i] /= a[i * d + i]; }       // L y = b
    for i in (0..d).rev() { for k in i + 1..d { z[i] -= a[k * d + i] * z[k]; } z[i] /= a[i * d + i]; } // Lᵀ z = y
    z
}

pub fn ridge(x: &[f64], y: &[f64], n: usize, d: usize, lam: f64) -> Vec<f64> {
    let a = xtwx(x, &vec![1.0; n], n, d, lam);
    let c: Vec<f64> = (0..d).map(|j| (0..n).map(|i| x[i * d + j] * y[i]).sum()).collect();
    cholesky_solve(a, &c, d)
}

fn sigmoid(z: f64) -> f64 { if z >= 0.0 { 1.0 / (1.0 + (-z).exp()) } else { let e = z.exp(); e / (1.0 + e) } }

/// Stable per-example log loss from the logit: softplus(z) − y·z.
fn log_loss(z: f64, y: f64) -> f64 { z.max(0.0) + (-z.abs()).exp().ln_1p() - y * z }

/// Logistic regression by Newton/IRLS. Returns (weights, iterations).
pub fn logistic_newton(x: &[f64], y: &[f64], n: usize, d: usize, lam: f64) -> (Vec<f64>, usize) {
    let mut w = vec![0.0; d];
    for it in 1..=50 {
        let z: Vec<f64> = (0..n).map(|i| (0..d).map(|j| x[i * d + j] * w[j]).sum()).collect();
        let p: Vec<f64> = z.iter().map(|&v| sigmoid(v)).collect();
        let wts: Vec<f64> = p.iter().map(|&q| (q * (1.0 - q)).max(1e-12)).collect();
        let mut g: Vec<f64> = (0..d).map(|j| (0..n).map(|i| x[i * d + j] * (p[i] - y[i])).sum()).collect();
        for j in 1..d { g[j] += lam * w[j]; }
        let step = cholesky_solve(xtwx(x, &wts, n, d, lam), &g, d);
        for j in 0..d { w[j] -= step[j]; }
        if step.iter().map(|s| s * s).sum::<f64>().sqrt() < 1e-10 { return (w, it); }
    }
    (w, 50)
}

fn main() {
    let (n, d) = (2000, 4);                                       // column 0 is the intercept
    let mut s = 0x2545F4914F6CDD1Du64;
    let mut u = || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s >> 11) as f64 / (1u64 << 53) as f64 };
    let true_w = [0.5, 2.0, -1.0, 0.0];
    let mut x = vec![0.0; n * d];
    let (mut yr, mut yc) = (vec![0.0; n], vec![0.0; n]);
    for i in 0..n {
        x[i * d] = 1.0;
        for j in 1..d { x[i * d + j] = 2.0 * u() - 1.0; }
        let z: f64 = (0..d).map(|j| x[i * d + j] * true_w[j]).sum();
        yr[i] = z + 0.1 * (u() - 0.5);                            // linear target + noise
        yc[i] = if u() < sigmoid(3.0 * z) { 1.0 } else { 0.0 };   // Bernoulli labels
    }
    println!("ridge w = {:.3?}", ridge(&x, &yr, n, d, 1e-3));
    let (w, iters) = logistic_newton(&x, &yc, n, d, 1.0);
    let loss: f64 = (0..n).map(|i| log_loss((0..d).map(|j| x[i * d + j] * w[j]).sum(), yc[i])).sum::<f64>() / n as f64;
    println!("logistic w = {:.3?} (truth 3x {:?}), {iters} Newton steps, log loss {loss:.4}", w, true_w);
}

Complexity and constants

• XᵀX: n·d²/2 multiply-adds (symmetric). With n = 10⁷, d = 100: 5·10¹⁰, order of a second on a multicore BLAS (syrk).
• Cholesky: d³/3 flops. d = 1000: order of 10-100 ms; d = 10⁴: order of seconds with a tuned multithreaded LAPACK.
• SGD/L-BFGS: O(n·d) per pass, memory-bound on X. Sparse features (hashed text): O(nnz) per pass. Vowpal Wabbit and liblinear train on order of 10⁸ sparse rows in minutes.
• Prediction: one dot product, about d multiply-adds.

Efficiency

• Cost: XᵀX is n·d²/2 multiply-adds (5·10¹⁰ for n = 10⁷, d = 100, order of a second with BLAS syrk); Cholesky adds d³/3 (d = 1000: order of 10-100 ms). Iterative solvers cost O(n·d) per pass, memory-bound on X; sparse data O(nnz). A prediction is d multiply-adds, tens of nanoseconds for d = 100 (order of).
• Levers:
  • Normal equations plus Cholesky — one pass over X, then solve in d³/3 — n ≫ d, moderate d, well-conditioned.
  • QR instead of Cholesky — conditioning κ instead of κ² (Golub and Van Loan) — ill-conditioned X, at about 2·n·d² cost.
  • Sparse storage plus feature hashing — cost per pass O(nnz), 10⁸ sparse rows in minutes (Vowpal Wabbit, liblinear; this page) — text and ad features.
  • Warm-started regularization path — a whole λ path for about the price of one fit (Friedman 2010) — tuning λ.
• Naive → standard → state of the art:
  • Invert XᵀX explicitly · O(n·d² + d³), numerically fragile
  • Cholesky or QR for least squares; L-BFGS or IRLS for logistic · one pass of X for the Gram matrix, tens of iterations for logistic (order of)
  • Sparse coordinate descent with warm starts, or online hashed training · O(nnz) per pass (this page)

Implementation notes

• Standardize features (zero mean, unit variance) before penalizing. Otherwise λ penalizes features by their units.
• Accumulate XᵀX in f64. In f32, sums over millions of rows lose precision and Cholesky can fail on a matrix that is mathematically positive definite.
• For ill-conditioned problems without regularization, solve the least-squares problem with QR (R w = Qᵀy) or SVD. Same answer, condition number κ instead of κ².
• Newton with a step-halving line search is robust against the rare non-decreasing step. Pure Newton (above) is fine with λ > 0.
• Class imbalance: use per-example weights in the gradient and Hessian, or fix the intercept afterwards by log(π_true/π_sample) when you subsample negatives.
• Sparse high-dimensional data: L1 (lasso) by coordinate descent with soft thresholding gives exact zeros (glmnet).

Pitfalls

• Perfect separation: unregularized logistic weights grow without bound and the solver reports non-convergence. Add L2.
• Multicollinearity: coefficients swing wildly and change sign between samples although predictions are fine. Do not interpret them. Use ridge.
• Penalizing the intercept: predictions shift toward 0.5 or toward 0.
• log(sigmoid(z)) computed literally as log(1/(1+exp(−z))): −inf once z < about −710 in f64 (−89 in f32). Use the softplus form.
• Reading coefficients as causal effects.
• Data leakage through standardization statistics computed on the full dataset including the test split.

Variants and alternatives

In production

scikit-learn (LinearRegression, Ridge, LogisticRegression via lbfgs, liblinear, saga), glmnet, statsmodels (IRLS with standard errors), liblinear, Vowpal Wabbit, Spark MLlib, ad click-through models (logistic regression on hashed feature crosses), statistical analysis in medicine and the social sciences (statsmodels, glmnet).

Implementations

Sources

• Hastie, Tibshirani, Friedman, The Elements of Statistical Learning, 2nd ed., Springer 2009, ch. 3-4. book
• Golub and Van Loan, Matrix Computations, 4th ed., 2013.
• Fan, Chang, Hsieh, Wang, Lin, LIBLINEAR: A Library for Large Linear Classification, JMLR 2008. paper · code
• Friedman, Hastie, Tibshirani, Regularization Paths for Generalized Linear Models via Coordinate Descent, Journal of Statistical Software 2010. paper
• McMahan et al., Ad Click Prediction: a View from the Trenches, KDD 2013. paper