From Scratch · C++ · ML · Regression · Linear Regression · Part 1

Ordinary Least Squares in C++

July 2026 · 30 min read

This post builds a reusable ordinary least squares (OLS) regressor in standard-library C++: no external linear-algebra dependency, only portable standard headers and vector<vector<double>>. It follows the data and model contract from the series overview: the caller supplies a raw feature matrix, the model owns the intercept, and every variant in this sub-series exposes the same fit / predict / evaluate lifecycle. The finished program fits a two-feature synthetic dataset and recovers the true coefficients to within a few hundredths, reaching roughly $0.98$ test $R^2$. The theory of the projection geometry behind this estimator is covered separately in OLS, Normal Equations, and Projection Geometry; this post is about turning that derivation into working, validated code.

The Least-Squares Objective

The training data is a feature matrix $X \in \mathbb{R}^{n \times d}$ ($n$ rows, one per sample; $d$ columns, one per feature) and a target vector $y \in \mathbb{R}^n$. A linear model predicts

$$\hat y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_d x_{id},$$

where $\beta_0$ is the intercept and $\beta_1, \dots, \beta_d$ are the feature coefficients. Ordinary least squares picks the $\beta$ that minimizes the sum of squared residuals:

$$J(\beta) = \sum_{i=1}^{n} (y_i - \hat y_i)^2 = \lVert y - X_b\beta \rVert_2^2,$$

where $X_b$ is the $n \times (d+1)$ design matrix obtained by prepending a column of ones to $X$, and $\beta \in \mathbb{R}^{d+1}$ stacks the intercept on top of the feature coefficients. Minimizing a sum of squares is a convex quadratic problem in $\beta$, so a single stationary point is also the global minimum whenever $X_b^\top X_b$ is invertible.

The Intercept Column

The series overview settles this once: the public interface never asks a caller for a bias column. OLSRegressor::fit receives the raw $n \times d$ matrix X and treats the intercept as an implicit column of ones prepended to it. That column contributes a fixed row/column of known values to $X_b^\top X_b$ and $X_b^\top y$, so it can be folded into the normal-equation assembly directly, without ever materializing $X_b$ in memory. That is the difference between an $n\times(d+1)$ allocation and reusing the caller's own $X$.

Deriving the Normal Equations

Setting the gradient of $J(\beta)$ to zero gives the normal equations. Expanding the squared norm,

$$J(\beta) = y^\top y - 2\beta^\top X_b^\top y + \beta^\top X_b^\top X_b \beta,$$

and differentiating with respect to $\beta$,

$$\nabla_\beta J(\beta) = -2X_b^\top y + 2X_b^\top X_b \beta = 0 \quad\Longrightarrow\quad X_b^\top X_b\,\beta = X_b^\top y.$$

Write $A = X_b^\top X_b$ (size $(d+1)\times(d+1)$, symmetric) and $b = X_b^\top y$ (size $d+1$). Solving $A\beta = b$ for $\beta$ is the entire fitting problem. $A$ is symmetric positive semi-definite always, and strictly positive definite (hence invertible) exactly when the columns of $X_b$ are linearly independent, which requires $n \ge d+1$ and no exact collinearity among the intercept column and the features.

Symbol	Shape	Meaning
$X$	$n \times d$	raw feature matrix, no bias column
$y$	$n$	regression targets
$X_b$	$n \times (d{+}1)$	$X$ with an implicit leading column of ones
$A = X_b^\top X_b$	$(d{+}1) \times (d{+}1)$	normal-equation matrix, symmetric
$b = X_b^\top y$	$d{+}1$	normal-equation right-hand side
$\beta$	$d{+}1$	$\beta_0$ intercept, $\beta_{1..d}$ coefficients

Avoiding Explicit Matrix Inversion

$\beta = A^{-1}b$ is the textbook formula, but computing $A^{-1}$ explicitly is the wrong move in code. Forming the full inverse costs about as much arithmetic as solving the system directly, produces a $(d+1)\times(d+1)$ matrix that is thrown away after one multiplication, and accumulates more rounding error than solving $A\beta = b$ directly, because every entry of $A^{-1}$ mixes information from the whole matrix instead of following one clean elimination pass. The implementation below solves the linear system in place with Gaussian elimination and never computes $A^{-1}$.

Condition number squares. A least-squares problem inherits its numerical sensitivity from $X_b$, but the normal equations work with $A = X_b^\top X_b$, and $\operatorname{cond}(A) = \operatorname{cond}(X_b)^2$. A design matrix with condition number $10^4$ (already uncomfortable) turns into a normal-equation matrix with condition number $10^8$. For well-scaled, non-collinear features this is fine; for near-collinear or poorly scaled features it can wipe out several digits of precision before elimination even starts. The next post replaces this step with QR and SVD based solvers that work directly on $X_b$ and never square the condition number.

Forming $A$ and $b$ Without Materializing $X_b$

Because the intercept column is always ones, every entry that touches it has a closed form: $A_{00} = n$, $A_{0,j+1} = A_{j+1,0} = \sum_i X_{ij}$, and $b_0 = \sum_i y_i$. The remaining feature-feature block is the usual Gram matrix, and $A$ is symmetric, so only its upper triangle needs computing. Matrix below is a plain vector<vector<double>> alias used throughout this implementation:

int p = d + 1;
Matrix A(p, vector<double>(p, 0.0));
vector<double> rhs(p, 0.0);

A[0][0] = (double)n;
for (int i = 0; i < n; ++i) rhs[0] += y[i];

for (int j = 0; j < d; ++j) {
    double colSum = 0.0, colDotY = 0.0;
    for (int i = 0; i < n; ++i) { colSum += X[i][j]; colDotY += X[i][j] * y[i]; }
    A[0][j + 1] += colSum;
    A[j + 1][0] += colSum;
    rhs[j + 1] += colDotY;
}
for (int j1 = 0; j1 < d; ++j1) {
    for (int j2 = j1; j2 < d; ++j2) {
        double s = 0.0;
        for (int i = 0; i < n; ++i) s += X[i][j1] * X[i][j2];
        A[j1 + 1][j2 + 1] += s;
        if (j2 != j1) A[j2 + 1][j1 + 1] += s;
    }
}

This costs $O(nd^2)$ arithmetic and $O(d^2)$ extra memory: no $n\times(d+1)$ design matrix is ever allocated, only the small $A$ and $b$.

A Partial-Pivoting Gaussian Solver

Solving $A\beta = b$ uses Gaussian elimination with partial pivoting: at each column, the row with the largest-magnitude entry in that column (among the remaining rows) is swapped into the pivot position before eliminating. Pivoting on the largest available magnitude keeps the elimination multipliers bounded and avoids dividing by a small number that would blow up rounding error.

vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0)
        throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size())
        throw invalid_argument("solveLinearSystem: A and b dimensions do not match");

    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " +
                                   to_string(i) + " has " + to_string(A[i].size()) +
                                   " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v))
                throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i]))
            throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() *
                       max(1.0, matrixScale) * n * 100.0;

    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i) {
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        }
        if (best <= tol) {
            throw runtime_error(
                "solveLinearSystem: matrix is singular or near-singular "
                "(no usable pivot in column " + to_string(k) + "); "
                "check for duplicated/collinear features");
        }
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }

        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }

    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

The function takes A and b by value on purpose: elimination mutates both in place, and the caller's normal-equation matrix should not be silently destroyed as a side effect of solving it. Because $(d+1)$ is small in most regression problems, this $O(d^3)$ elimination plus $O(d^2)$ back-substitution is cheap relative to the $O(nd^2)$ cost of forming $A$ and $b$ in the first place.

The walkthrough below steps through elimination on a tiny three-parameter system (intercept plus two features) built from five data points, so the pivot search, the row swap, and the two elimination passes are all visible before back-substitution recovers $\beta$.

▶ Solving $A\beta = b$ by Partial-Pivoting Elimination

Columns are $\beta_0$ (intercept), $\beta_1$, $\beta_2$, and the right-hand side $b$. Rows are the three normal equations. Each step shows the pivot search selecting the largest-magnitude entry in a column before elimination below it.

Validation and Errors

fit checks its inputs before touching the normal equations, and every failure is an explicit, specific exception rather than a silently wrong answer:

Row count mismatch between X and y, or a ragged X (rows with different lengths), throws std::invalid_argument naming the offending row.
An empty training set or zero feature columns throws std::invalid_argument.
Non-finite entries (NaN or inf) in X or y throw std::invalid_argument naming the row, since a single bad value would otherwise silently poison every sum in $A$ and $b$.
Fewer than $d+1$ rows (not enough data to determine $d+1$ parameters) throws std::invalid_argument before the solver is even invoked.
The standalone solver rejects an empty or non-square $A$, a mismatched $b$, and non-finite matrix entries with std::invalid_argument.
A singular or near-singular $A$ (perfectly collinear features, or a constant column that duplicates the intercept) throws std::runtime_error from inside solveLinearSystem.
Calling predict on an unfitted model throws std::logic_error; calling it with the wrong feature count throws std::invalid_argument.
The metric functions reject mismatched, empty, or non-finite vectors; r2 also rejects a constant target because $R^2$ is undefined when target variance is zero.

None of these are caught and papered over inside the library code. A caller that wants a fallback (dropping a redundant column, imputing a missing value) makes that decision explicitly at its own call site; the regressor's job is only to refuse to silently fit garbage.

Fit and Predict

Putting the pieces together, OLSRegressor validates its inputs, assembles $A$ and $b$, solves for $\beta$, and splits the result back into an intercept and a coefficient vector:

struct OLSRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;

    void fit(const Matrix& X, const vector<double>& y) {
        if (X.size() != y.size())
            throw invalid_argument("OLSRegressor::fit: X has " + to_string(X.size()) +
                                    " rows but y has " + to_string(y.size()) + " entries");
        int n = (int)X.size();
        if (n == 0) throw invalid_argument("OLSRegressor::fit: empty training set");

        int d = (int)X[0].size();
        if (d == 0) throw invalid_argument("OLSRegressor::fit: zero feature columns");
        for (int i = 0; i < n; ++i) {
            if ((int)X[i].size() != d)
                throw invalid_argument("OLSRegressor::fit: row " + to_string(i) +
                                        " has " + to_string(X[i].size()) +
                                        " columns, expected " + to_string(d));
            for (double v : X[i])
                if (!isfinite(v))
                    throw invalid_argument("OLSRegressor::fit: non-finite feature value in row " + to_string(i));
            if (!isfinite(y[i]))
                throw invalid_argument("OLSRegressor::fit: non-finite target at row " + to_string(i));
        }
        if (n < d + 1)
            throw invalid_argument("OLSRegressor::fit: need at least " + to_string(d + 1) +
                                    " rows for " + to_string(d) + " features plus intercept, got " +
                                    to_string(n));

        int p = d + 1;
        Matrix A(p, vector<double>(p, 0.0));
        vector<double> rhs(p, 0.0);

        A[0][0] = (double)n;
        for (int i = 0; i < n; ++i) rhs[0] += y[i];

        for (int j = 0; j < d; ++j) {
            double colSum = 0.0, colDotY = 0.0;
            for (int i = 0; i < n; ++i) { colSum += X[i][j]; colDotY += X[i][j] * y[i]; }
            A[0][j + 1] += colSum;
            A[j + 1][0] += colSum;
            rhs[j + 1] += colDotY;
        }
        for (int j1 = 0; j1 < d; ++j1) {
            for (int j2 = j1; j2 < d; ++j2) {
                double s = 0.0;
                for (int i = 0; i < n; ++i) s += X[i][j1] * X[i][j2];
                A[j1 + 1][j2 + 1] += s;
                if (j2 != j1) A[j2 + 1][j1 + 1] += s;
            }
        }

        vector<double> beta = solveLinearSystem(A, rhs);
        intercept = beta[0];
        coef.assign(beta.begin() + 1, beta.end());
        nFeatures = d;
        fitted = true;
    }

    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("OLSRegressor::predict: model is not fitted yet");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("OLSRegressor::predict: expected " + to_string(nFeatures) +
                                    " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) out += coef[j] * x[j];
        return out;
    }

    vector<double> predict(const Matrix& X) const {
        vector<double> out;
        out.reserve(X.size());
        for (const auto& row : X) out.push_back(predict(row));
        return out;
    }
};

Evaluation Metrics

The same three free functions will be reused by every later post in the sub-series, so they live outside the regressor:

void validateMetricInputs(const vector<double>& yTrue,
                          const vector<double>& yPred,
                          const string& name) {
    if (yTrue.size() != yPred.size())
        throw invalid_argument(name + ": size mismatch");
    if (yTrue.empty())
        throw invalid_argument(name + ": vectors must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i) {
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(name + ": non-finite value at index " + to_string(i));
    }
}

double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "mse");
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        double e = yTrue[i] - yPred[i];
        s += e * e;
    }
    return s / yTrue.size();
}

double rmse(const vector<double>& yTrue, const vector<double>& yPred) {
    return sqrt(mse(yTrue, yPred));
}

double r2(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "r2");
    double meanY = accumulate(yTrue.begin(), yTrue.end(), 0.0) / yTrue.size();
    double ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i] - yPred[i]) * (yTrue[i] - yPred[i]);
        ssTot += (yTrue[i] - meanY) * (yTrue[i] - meanY);
    }
    if (ssTot == 0.0)
        throw invalid_argument("r2: target values are constant, so R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}

rmse reports error in the same units as $y$; $R^2$ compares residual error with a mean-only baseline, with $1.0$ a perfect fit, $0.0$ matching that baseline, and negative values possible for a worse fit. A constant target has no variance, so $R^2$ is undefined and the function throws instead of dividing by zero.

Dry Run on a Synthetic Dataset

The complete program below generates $240$ points from $y = 3 + 2x_1 - 1.5x_2 + \varepsilon$ with Gaussian noise, splits them 80/20, fits OLSRegressor on the training rows, and reports the recovered parameters plus train/test metrics. It also exercises the two error paths described above so the exception messages can be seen directly instead of only described:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;

using Matrix = vector<vector<double>>;

vector<double> solveLinearSystem(Matrix A, vector<double> b) {
    int n = (int)A.size();
    if (n == 0)
        throw invalid_argument("solveLinearSystem: A must be non-empty");
    if (b.size() != A.size())
        throw invalid_argument("solveLinearSystem: A and b dimensions do not match");

    double matrixScale = 0.0;
    for (int i = 0; i < n; ++i) {
        if (A[i].size() != A.size())
            throw invalid_argument("solveLinearSystem: A must be square; row " +
                                   to_string(i) + " has " + to_string(A[i].size()) +
                                   " entries, expected " + to_string(n));
        for (double v : A[i]) {
            if (!isfinite(v))
                throw invalid_argument("solveLinearSystem: A contains a non-finite value");
            matrixScale = max(matrixScale, fabs(v));
        }
        if (!isfinite(b[i]))
            throw invalid_argument("solveLinearSystem: b contains a non-finite value");
    }
    const double tol = numeric_limits<double>::epsilon() *
                       max(1.0, matrixScale) * n * 100.0;

    for (int k = 0; k < n; ++k) {
        int piv = k;
        double best = fabs(A[k][k]);
        for (int i = k + 1; i < n; ++i) {
            if (fabs(A[i][k]) > best) { best = fabs(A[i][k]); piv = i; }
        }
        if (best <= tol) {
            throw runtime_error(
                "solveLinearSystem: matrix is singular or near-singular "
                "(no usable pivot in column " + to_string(k) + "); "
                "check for duplicated/collinear features");
        }
        if (piv != k) { swap(A[k], A[piv]); swap(b[k], b[piv]); }

        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            if (factor == 0.0) continue;
            for (int j = k; j < n; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }

    vector<double> beta(n, 0.0);
    for (int k = n - 1; k >= 0; --k) {
        double s = b[k];
        for (int j = k + 1; j < n; ++j) s -= A[k][j] * beta[j];
        beta[k] = s / A[k][k];
    }
    return beta;
}

struct OLSRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;

    void fit(const Matrix& X, const vector<double>& y) {
        if (X.size() != y.size())
            throw invalid_argument("OLSRegressor::fit: X has " + to_string(X.size()) +
                                    " rows but y has " + to_string(y.size()) + " entries");
        int n = (int)X.size();
        if (n == 0) throw invalid_argument("OLSRegressor::fit: empty training set");

        int d = (int)X[0].size();
        if (d == 0) throw invalid_argument("OLSRegressor::fit: zero feature columns");
        for (int i = 0; i < n; ++i) {
            if ((int)X[i].size() != d)
                throw invalid_argument("OLSRegressor::fit: row " + to_string(i) +
                                        " has " + to_string(X[i].size()) +
                                        " columns, expected " + to_string(d));
            for (double v : X[i])
                if (!isfinite(v))
                    throw invalid_argument("OLSRegressor::fit: non-finite feature value in row " + to_string(i));
            if (!isfinite(y[i]))
                throw invalid_argument("OLSRegressor::fit: non-finite target at row " + to_string(i));
        }
        if (n < d + 1)
            throw invalid_argument("OLSRegressor::fit: need at least " + to_string(d + 1) +
                                    " rows for " + to_string(d) + " features plus intercept, got " +
                                    to_string(n));

        int p = d + 1;
        Matrix A(p, vector<double>(p, 0.0));
        vector<double> rhs(p, 0.0);

        A[0][0] = (double)n;
        for (int i = 0; i < n; ++i) rhs[0] += y[i];

        for (int j = 0; j < d; ++j) {
            double colSum = 0.0, colDotY = 0.0;
            for (int i = 0; i < n; ++i) { colSum += X[i][j]; colDotY += X[i][j] * y[i]; }
            A[0][j + 1] += colSum;
            A[j + 1][0] += colSum;
            rhs[j + 1] += colDotY;
        }
        for (int j1 = 0; j1 < d; ++j1) {
            for (int j2 = j1; j2 < d; ++j2) {
                double s = 0.0;
                for (int i = 0; i < n; ++i) s += X[i][j1] * X[i][j2];
                A[j1 + 1][j2 + 1] += s;
                if (j2 != j1) A[j2 + 1][j1 + 1] += s;
            }
        }

        vector<double> beta = solveLinearSystem(A, rhs);
        intercept = beta[0];
        coef.assign(beta.begin() + 1, beta.end());
        nFeatures = d;
        fitted = true;
    }

    double predict(const vector<double>& x) const {
        if (!fitted) throw logic_error("OLSRegressor::predict: model is not fitted yet");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("OLSRegressor::predict: expected " + to_string(nFeatures) +
                                    " features, got " + to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) out += coef[j] * x[j];
        return out;
    }

    vector<double> predict(const Matrix& X) const {
        vector<double> out;
        out.reserve(X.size());
        for (const auto& row : X) out.push_back(predict(row));
        return out;
    }
};

void validateMetricInputs(const vector<double>& yTrue,
                          const vector<double>& yPred,
                          const string& name) {
    if (yTrue.size() != yPred.size())
        throw invalid_argument(name + ": size mismatch");
    if (yTrue.empty())
        throw invalid_argument(name + ": vectors must be non-empty");
    for (size_t i = 0; i < yTrue.size(); ++i) {
        if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
            throw invalid_argument(name + ": non-finite value at index " + to_string(i));
    }
}

double mse(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "mse");
    double s = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        double e = yTrue[i] - yPred[i];
        s += e * e;
    }
    return s / yTrue.size();
}

double rmse(const vector<double>& yTrue, const vector<double>& yPred) {
    return sqrt(mse(yTrue, yPred));
}

double r2(const vector<double>& yTrue, const vector<double>& yPred) {
    validateMetricInputs(yTrue, yPred, "r2");
    double meanY = accumulate(yTrue.begin(), yTrue.end(), 0.0) / yTrue.size();
    double ssRes = 0.0, ssTot = 0.0;
    for (size_t i = 0; i < yTrue.size(); ++i) {
        ssRes += (yTrue[i] - yPred[i]) * (yTrue[i] - yPred[i]);
        ssTot += (yTrue[i] - meanY) * (yTrue[i] - meanY);
    }
    if (ssTot == 0.0)
        throw invalid_argument("r2: target values are constant, so R^2 is undefined");
    return 1.0 - ssRes / ssTot;
}

int main() {
    mt19937 rng(7);
    uniform_real_distribution<double> ux(-3.0, 3.0);
    normal_distribution<double> noise(0.0, 0.6);

    const int N = 240;
    const double trueIntercept = 3.0, trueB1 = 2.0, trueB2 = -1.5;

    Matrix X;
    vector<double> y;
    for (int i = 0; i < N; ++i) {
        double x1 = ux(rng), x2 = ux(rng);
        X.push_back({x1, x2});
        y.push_back(trueIntercept + trueB1 * x1 + trueB2 * x2 + noise(rng));
    }

    vector<int> order(N);
    iota(order.begin(), order.end(), 0);
    shuffle(order.begin(), order.end(), rng);

    Matrix Xtr, Xte;
    vector<double> ytr, yte;
    int nTrain = (int)(0.8 * N);
    for (int k = 0; k < N; ++k) {
        int i = order[k];
        if (k < nTrain) { Xtr.push_back(X[i]); ytr.push_back(y[i]); }
        else            { Xte.push_back(X[i]); yte.push_back(y[i]); }
    }

    OLSRegressor model;
    model.fit(Xtr, ytr);

    printf("intercept = %.4f (true %.1f)\n", model.intercept, trueIntercept);
    printf("coef[0]   = %.4f (true %.1f)\n", model.coef[0], trueB1);
    printf("coef[1]   = %.4f (true %.1f)\n", model.coef[1], trueB2);

    auto predTr = model.predict(Xtr);
    auto predTe = model.predict(Xte);
    printf("train RMSE = %.4f   train R^2 = %.4f\n", rmse(ytr, predTr), r2(ytr, predTr));
    printf("test  RMSE = %.4f   test  R^2 = %.4f\n", rmse(yte, predTe), r2(yte, predTe));

    // Error-handling demo: mismatched dimensions.
    try {
        Matrix badX = {{1.0, 2.0}, {3.0}};
        OLSRegressor bad;
        bad.fit(badX, {1.0, 2.0});
    } catch (const invalid_argument& e) {
        printf("caught expected error: %s\n", e.what());
    }

    // Error-handling demo: singular normal equations (duplicated column).
    try {
        Matrix dupX = {{1.0, 1.0}, {2.0, 2.0}, {3.0, 3.0}};
        vector<double> dupY = {1.0, 2.0, 3.0};
        OLSRegressor dup;
        dup.fit(dupX, dupY);
    } catch (const runtime_error& e) {
        printf("caught expected error: %s\n", e.what());
    }
    return 0;
}

Compiled with g++ -O2 -std=c++17, this prints:

intercept = 2.9639 (true 3.0) coef[0] = 1.9777 (true 2.0) coef[1] = -1.5173 (true -1.5) train RMSE = 0.6111 train R^2 = 0.9774 test RMSE = 0.5638 test R^2 = 0.9813 caught expected error: OLSRegressor::fit: row 1 has 1 columns, expected 2 caught expected error: solveLinearSystem: matrix is singular or near-singular (no usable pivot in column 2); check for duplicated/collinear features

The recovered intercept and coefficients land close to the true generating values $3.0$, $2.0$, and $-1.5$; the gap is exactly what the injected noise ($\sigma = 0.6$) should leave behind. Test $R^2$ sits slightly above train $R^2$ here only because of how this particular random split happened to land, not because the model generalizes better than it fits; with more data both numbers converge to the same value.

Complexity and Memory

Forming $A$ and $b$ costs $O(nd^2)$ time while reusing $X$ in place; elimination and back-substitution cost $O(d^3)$ and $O(d^2)$ respectively. Memory beyond the caller's own $X$ and $y$ is $O(d^2)$, the size of $A$, plus $O(d)$ for $b$ and $\beta$. In particular, no $O(nd)$ design-matrix copy is ever allocated. When $n \gg d$, which is the common regression regime, the $O(nd^2)$ accumulation dominates the $O(d^3)$ solve; when $d$ approaches $n$, the solve cost catches up and, per the row-count check in fit, $n$ must stay at least $d+1$ for the system to be solvable at all.

Limitations and Edge Cases

Perfect collinearity. A duplicated or exactly linearly dependent feature column makes $A$ singular; solveLinearSystem throws rather than returning a coefficient vector built on an arbitrary pseudo-inverse choice.
Near collinearity. Strongly correlated features do not trip the singularity check but do inflate $\operatorname{cond}(A) = \operatorname{cond}(X_b)^2$, so small changes in the data can swing $\beta$ noticeably. This is exactly the situation the QR and SVD solvers post addresses without squaring the condition number.
Too few rows. $n \le d$ leaves the system underdetermined; fit refuses before attempting a solve.
Outliers. Squared loss weights large residuals quadratically, so a handful of extreme points can dominate the fit. The robust and quantile regression post replaces the loss rather than the solver to address this.
Constant targets. OLS can fit a constant target through the intercept, but $R^2$ is undefined because the target has zero variance; the metric function throws an explicit std::invalid_argument.
Non-linear relationships. This estimator is linear in the parameters; a curved relationship in $X$ needs transformed or additional features before OLS can capture it, which is a data question, not something the solver can fix.

Next Step

The normal equations are the right tool when $X_b$ is reasonably well-conditioned and $d$ is modest, since the whole fit reduces to one small linear solve. When the condition number of $X_b$ is already large, squaring it inside $A$ can cost real accuracy before elimination even begins. The next post keeps the same public interface, fit and predict with the identical shape contract, and swaps the normal-equation solve for QR and SVD based solvers that operate on $X_b$ directly.