Regression from Scratch in C++

July 2026 · Series overview

This gateway joins two model families that share a linear score but answer different questions. Linear regression estimates a continuous conditional mean. Logistic regression turns the same score into a Bernoulli probability and then, only when a decision is needed, into a class. Every implementation uses portable C++17 and the standard library only. Matrices remain vector<vector<double>>, the intercept remains model-owned, and numerical failures remain explicit rather than silently hidden.

One Linear Score, Two Statistical Models

Both branches begin with $z_i=\beta_0+x_i^\top\beta$. The interpretation after that point is different:

The word regression in logistic regression describes a linear model for log-odds. Its observed target is still categorical. Treating a binary label as an unrestricted continuous response can produce values below zero or above one; applying a sigmoid repairs the range, while the Bernoulli likelihood supplies the matching objective.

Shared Data and Model Contracts

using Matrix = std::vector<std::vector<double>>;

// X has n rows and d finite features per row. No bias column is supplied.
// The model owns beta_0 and stores d feature coefficients.
// fit validates all rows before changing fitted state.

All implementations preserve a strong lifecycle. Construction creates an unfitted object. fit validates shapes, finiteness, targets, and solver settings before committing parameters. Prediction before fitting raises std::logic_error. Invalid caller data raises std::invalid_argument. A numerical solver that cannot make a reliable step raises std::runtime_error. Failed fitting leaves no partly fitted model behind.

The classifier adds three distinct operations: predictLogit returns $z$, predictProba returns $\sigma(z)$, and predict applies a validated threshold. Keeping these operations separate prevents a threshold policy from being confused with probability estimation.

Intercept and Train-Only Scaling

The caller never appends a column of ones. Each solver accounts for the intercept internally, so preprocessing cannot accidentally scale it or add it twice. Feature scaling is estimated on training rows only and then frozen. The same constants transform validation, test, and production rows. This rule avoids leakage and improves the geometry seen by gradient descent and Newton methods.

A constant feature needs explicit treatment because its scale is zero. The linear-series gradient implementation maps its scale to one after centering, leaving a zero column. Logistic implementations may reject an unidentifiable design or rely on regularization, but no implementation divides by zero or silently produces a non-finite parameter.

Numerical Reliability as Part of the API

• Stable primitives: logistic loss uses softplus from logits, not log(sigmoid(z)); the sigmoid branches on the sign of $z$.
• No explicit inverse: direct solvers solve linear systems with pivoting. QR or SVD is preferred when normal equations are poorly conditioned.
• Monotonic safeguards: iterative solvers track finite loss, reject non-finite steps, and use convergence tolerances tied to objective or parameter change.
• Observable history: iteration records expose progress and make stalled or divergent training testable.
• Finite boundaries: every feature, target, parameter, threshold, and reported metric is checked before it crosses the public boundary.

A Target-to-Diagnostics Decision Workflow

The workflow begins with the target support, not with a favorite solver. Each step reveals the estimand, matching loss, practical solver family, and diagnostics for both branches. The status line also fits the stated six-row examples: closed-form simple OLS supplies the continuous RMSE, while deterministic batch gradient descent supplies the Bernoulli log-loss. Those computed values make the diagram an executable model-selection map rather than a decorative flowchart.

Published Roadmap

Implementation and Theory Map

The implementation curriculum and the mathematical curriculum reinforce one another. The linear branch pairs its C++ pages with the theory sequence on OLS and projection, numerical solvers, gradient descent, regularization, robust and quantile methods, and diagnostics. The logistic branch pairs stable primitives and both optimizers with likelihood and log-loss and gradient, Newton, and IRLS theory; later implementation modules align with the theory pages on regularization and separation, softmax, classification metrics, and calibration.

Code reuse follows the same boundaries. Matrix and validation helpers, train-only scalers, finite checks, deterministic data, metric accumulators, and linear-system solvers remain reusable. The response model remains specific: squared loss and residual analysis belong to continuous targets, while sigmoid probabilities, Bernoulli loss, calibration, and threshold policy belong to binary targets.

Choosing a Route

A modest, well-conditioned continuous problem can begin with OLS. Correlated columns motivate QR, SVD, or regularization. Large row counts motivate gradient methods. Binary targets begin with stable log-loss; a first-order gradient method offers simple $O(nd)$ iterations, while Newton/IRLS spends $O(nd^2+d^3)$ per iteration to use curvature and often needs far fewer iterations. Perfectly separated classes require regularization or an explicit bounded stopping policy because the unpenalized finite optimum does not exist.

Verification Across the Series

Each complete listing is compiled with strict C++17 warnings and exercised on deterministic data. Property checks complement printed examples: probabilities stay in $[0,1]$, stable loss remains finite for extreme logits, analytic gradients agree with central differences, accepted steps do not increase the objective, and predictions fail clearly when lifecycle or shape contracts are violated. These checks make the derivations executable rather than decorative.