Linear Regression from Scratch in C++

July 2026 · Series overview

This sub-series builds the linear-regression family in standard-library C++: no external linear-algebra dependency, only portable standard headers and vector<vector<double>>. It starts from the closed-form ordinary least squares solution and extends outward to numerically robust solvers, iterative optimization, regularized variants, robust losses, and diagnostics. The companion Machine Learning · Linear Regression series covers the conceptual treatment, projection geometry, statistical assumptions, and inference; this sub-series turns that theory into working, validated code.

Series Scope

Every post here compiles with a standard C++17 compiler and nothing else. There is no BLAS, no Eigen, no external matrix library: matrices are vector<vector<double>>, and every numerical routine, elimination, decomposition, or gradient step, is written out explicitly. That is a deliberate trade: production code should reach for a tuned linear-algebra library, but writing the routines by hand is how the normal equations, QR, SVD, and gradient-based solvers stop being black boxes.

The roadmap below covers six posts: the closed-form solution, two numerically safer solvers, an iterative alternative that scales to large $n$, penalized variants for collinear or high-dimensional data, robust and quantile variants for non-Gaussian errors, and a shared set of diagnostics that apply across all of them.

The Data and Model Contract

Every post in this sub-series fits the same shape of problem, so the interface never changes even as the solver underneath does:

using Matrix = vector<vector<double>>;

// X: n rows (samples) x d columns (features). No bias column.
// y: n regression targets, one per row of X.
// Every row of X must have exactly d entries, and n must equal y.size().

A prediction is $\hat y = \beta_0 + \sum_{j=1}^{d}\beta_j x_j$ for every variant in the series; what changes from post to post is how $\beta$ is estimated, not what it means once found.

Handling the Intercept

The intercept is handled the same way in every post: the model owns it. Callers pass a raw $n \times d$ matrix with no column of ones; each fit implementation augments the problem internally, either by folding an implicit ones-column into the normal equations directly, as the ordinary least squares post does, or by centering $X$ and $y$ before optimizing and recovering the intercept from the means afterward, as the gradient-based and regularized posts do. Either way, the public contract is identical: fit(X, y) takes only the raw feature matrix, and the fitted model exposes a separate intercept alongside a coef vector of length $d$.

Keeping the intercept out of the public matrix avoids a common bug class: an intercept column of ones that a caller forgets to add, adds twice, or scales along with the other features during preprocessing. Because the series' data-preprocessing sub-series scales and encodes features before they reach any model here, an implicit intercept means preprocessing never has to special-case a bias column.

A Coherent Standard-Library Architecture

Every regressor in the series is a plain struct with two public entry points and a small set of validated fields, matching the pattern already used for the regression tree and GBM posts elsewhere in the from-scratch series:

struct SomeRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;          // length nFeatures once fitted

    void fit(const Matrix& X, const vector<double>& y);
    double predict(const vector<double>& x) const;
    vector<double> predict(const Matrix& X) const;
};

Exceptions follow one convention across every post: std::invalid_argument for a dimension or shape problem the caller could have caught before calling, std::runtime_error for a solver that failed on otherwise well-shaped input (a singular matrix, a non-converging iteration), and std::logic_error for calling predict before fit. None of the posts in this series swallow an exception with a broad catch or silently fall back to a default; a caller that wants a fallback makes that decision explicitly at its own call site.

The fit / predict / evaluate Lifecycle

Every post follows the same three-stage lifecycle:

• fit(X, y) validates the shapes, estimates $\beta_0$ and $\beta_{1..d}$ by whatever method that post covers, closed-form, iterative, or penalized, and marks the model fitted.
• predict(x) or the batch overload predict(X) applies the fitted linear form to new rows, after checking the model is fitted and the feature count matches.
• evaluate is not a member function; it is a small set of free functions, mse, rmse, and r2, shared across every post in the series, plus the residual and leverage diagnostics added in the final post.

Because every post shares this lifecycle, swapping OLSRegressor for a ridge or gradient-descent regressor later in the series is a one-line change at the call site: the shape of X, the shape of y, and the meaning of intercept and coef never move.

Solver and Variant Roadmap

Numerical Reliability

The normal equations square the condition number: $\operatorname{cond}(X_b^\top X_b) = \operatorname{cond}(X_b)^2$. For well-scaled features with no severe collinearity and a modest number of parameters, that is a non-issue and the direct solve in the OLS post is both the simplest and the fastest option. As features become more correlated, more poorly scaled, or more numerous relative to the number of rows, the squared condition number can erase real precision before elimination even starts. That is the entire motivation for the QR and SVD solvers post: same public interface, same $\hat y = \beta_0 + \sum_j \beta_j x_j$ prediction rule, a numerically steadier way of finding $\beta$.

Why Each Later Post Matters