QR and SVD Least-Squares Solvers in C++

July 2026 · 55 min read

The previous post solved least squares by assembling the $(d{+}1) \times (d{+}1)$ normal-equation matrix $A = X_b^\top X_b$ and applying Gaussian elimination to $A\beta = X_b^\top y$. That is the right approach for well-scaled, non-collinear features, but it pays a precision tax that grows quadratically with the condition of $X_b$. This post keeps the same public interface, the same fit / predict / evaluate contract, and the same fit on raw $X$ with an implicit intercept, and replaces the normal-equation solve with two solvers that operate directly on the augmented design matrix $X_b$: a Householder QR decomposition that handles full-column-rank designs, and a one-sided Jacobi SVD that additionally produces a minimum-norm solution for rank-deficient designs. All code uses only standard-library headers and vector<vector<double>>.

Conditioning and the Normal Equations

Every finite-precision computation has a condition number story. For least squares the relevant quantity is $\operatorname{cond}(X_b)$, the ratio of the largest to smallest singular value of the augmented design matrix. A condition number of $10^4$ means a $10^4$-fold amplification of input perturbations in the worst case. The problem is that the normal equations work with $A = X_b^\top X_b$, and forming this product squares the condition number:

A moderately ill-conditioned $X_b$ with $\operatorname{cond}(X_b) = 10^4$ produces an $A$ with $\operatorname{cond}(A) = 10^8$. Double-precision arithmetic carries about $15{-}16$ decimal digits, so $10^8$ conditioning leaves only $7{-}8$ reliable digits in the solution. For many regression problems this is acceptable; for near-collinear or poorly scaled features it can wipe out the signal before the solver even runs. The two solvers introduced here avoid this squaring entirely: QR and SVD operate on $X_b$ itself, so they only inherit $\operatorname{cond}(X_b)$, not its square.

The Augmented Design Matrix

Both solvers augment the raw $n \times d$ feature matrix $X$ with a leading column of ones to produce the $n \times (d{+}1)$ augmented matrix $X_b$, exactly as the normal-equation formulation does conceptually. Throughout this post $p = d{+}1$ denotes the number of parameters (intercept plus $d$ coefficients). The augmented matrix must be tall: $n \ge p$ is required for any unique least-squares solution to exist.

Materializing $X_b$ as an explicit $n \times p$ matrix costs $O(np)$ memory. Both solvers here allocate that matrix internally, unlike the normal-equation approach in the previous post, which avoided it entirely. This is the trade: the stability gain of working directly on $X_b$ comes at the cost of one additional $O(np)$ allocation on top of the caller's own $X$.

Householder QR

The QR decomposition factors $X_b = QR$ where $Q \in \mathbb{R}^{n \times n}$ is orthogonal ($Q^\top Q = I$) and $R \in \mathbb{R}^{n \times p}$ is upper triangular. For a tall full-rank matrix, only the first $p$ rows of $R$ are nonzero, so the least-squares problem reduces to:

where $R_p$ is the top $p \times p$ block of $R$ and $(Q^\top y)_{0:p}$ is the first $p$ entries of $Q^\top y$. Since $R_p$ is upper triangular and (for a full-rank design) nonsingular, the solution follows directly from back-substitution. No squaring of the condition number occurs because the orthogonal factor $Q$ preserves norms exactly.

Householder Reflectors

Householder triangularization applies a sequence of orthogonal reflections. At step $k$ the goal is to zero out all entries in column $k$ below the diagonal. The reflector is constructed as follows. Let $x = X_b[k:n, k]$ (the sub-column starting at the diagonal). Define

where $e_1$ is the first standard basis vector of the same length as $x$. The sign convention for $\alpha$ is chosen to avoid catastrophic cancellation when $x_0$ is close to $\lVert x \rVert_2$. The Householder matrix is

Multiplying $x$ by $H$ maps it to $\alpha e_1$: the entry at position $(k,k)$ becomes $\alpha$ (nonzero by construction) and all entries below it become zero. Because $H$ is orthogonal, $\lVert Hx \rVert_2 = \lVert x \rVert_2$, so no amplification of errors occurs. In code the reflector is never formed explicitly: the action of $H$ on any vector $u$ is computed as $Hu = u - \tau v$ where $\tau = 2(v^\top u)/(v^\top v)$, reducing the work from $O(n^2)$ for a dense matrix multiplication to $O(n)$ per application.

The algorithm applies the same reflector to all remaining columns of $X_b$ in column $k$ onward, and simultaneously to the running right-hand side $Q^\top y$ (initialized to $y$ and updated in place). After $p$ steps, $X_b$ has been overwritten with $R$ (upper triangular) and $Q^\top y$ holds the rotated right-hand side. Only the first $p$ entries of $Q^\top y$ are needed for the solve; the remaining $n{-}p$ entries are the norm of the residual but not the solution.

One Reflector by Hand

A single numeric pass removes any ambiguity from the reflector construction. The smallest instructive augmented matrix is $X_b = \begin{bmatrix}1&1\\1&2\\1&4\end{bmatrix}$ with intercept column $b_0 = [1,1,1]^\top$, feature column $b_1 = [1,2,4]^\top$, and right-hand side $y = [1,3,7]^\top$. The first reflector ($k=0$) targets the intercept sub-column $x = [1,1,1]^\top$ and drives its below-diagonal entries to zero. Every quantity below is exact to four decimals.

1. Sub-column norm. $\sigma = \lVert x \rVert_2 = \sqrt{1^2+1^2+1^2} = \sqrt{3} \approx 1.7321$. This is the magnitude the reflected pivot entry will carry, because an orthogonal reflection preserves length.
2. Stable sign. The pivot entry is $x_0 = 1 > 0$, so $\alpha = -\operatorname{sign}(x_0)\,\sigma = -1.7321$. Choosing the sign opposite to $x_0$ makes the leading reflector entry $v_0 = x_0 - \alpha$ a sum of like-signed terms rather than a difference, which is the cancellation-avoiding choice.
3. Reflector vector and its squared norm. $v = x - \alpha e_1 = [\,1-(-1.7321),\,1,\,1\,]^\top = [2.7321,\,1,\,1]^\top$, and $v^\top v = 2.7321^2 + 1^2 + 1^2 = 9.4641$. The reflector is applied through $v$ and $v^\top v$ alone; the dense matrix $H$ is never formed.
4. Projection scalar and reflected pivot column. Applying $H u = u - \tau v$ with $\tau = 2\,(v^\top u)/(v^\top v)$ to the pivot column itself gives $v^\top x = 4.7321$, so $\tau = 2(4.7321)/9.4641 = 1.0000$ and $x - \tau v = [-1.7321,\,0,\,0]^\top = \alpha e_1$. The pivot entry becomes exactly $\alpha$ and both entries below it vanish.
5. Same reflector on the feature column. For $b_1 = [1,2,4]^\top$ the projection scalar is $v^\top b_1 = 8.7321$, so $\tau = 1.8453$ and the reflected column is $[-4.0415,\,0.1547,\,2.1547]^\top$. The first entry fills the top row of $R$; the two lower entries become the sub-column that the next reflector will process.
6. Same reflector on the right-hand side. The identical reflector is applied to $y = [1,3,7]^\top$: $v^\top y = 12.7321$, $\tau = 2.6906$, and the rotated right-hand side is $[-6.3509,\,0.3094,\,4.3094]^\top$. This keeps $Q^\top y$ synchronized with $R$ without ever materializing $Q$.

A second reflector on the lower sub-column $[0.1547,\,2.1547]^\top$ finishes the triangularization, and back-substitution on the resulting $2 \times 2$ system recovers $\beta_1 = 2$ and $\beta_0 = -1$ exactly. The animation carries this same $3 \times 2$ example through both reflections and the back-substitution, one orthogonal rotation at a time.

Rank Checks Without Column Pivoting

Householder QR without column pivoting does not rearrange columns, so a rank-deficient design may not produce a clean zero on the diagonal: rounding can place a small but nonzero value there instead. After all $p$ reflections, the implementation scans the diagonal entries $|R_{kk}|$ for $k = 0, \dots, p{-}1$ and counts those above a tolerance:

If the count falls below $p$, the design is treated as rank-deficient and a std::runtime_error is thrown with a clear message directing the caller to the SVD solver. This is a deliberate design choice: QR without pivoting cannot reliably compute a minimum-norm solution for rank-deficient problems (the pivot column order matters for that), so it is better to fail explicitly and let the caller escalate to SVD.

For full-rank designs the condition estimate is the ratio of the largest to smallest diagonal magnitude:

This is a rough approximation of the true condition number of $X_b$ (not exact because the diagonal entries of $R$ are not the singular values), but it is cheap to compute and captures large differences in column scaling.

QR Regressor Implementation

QRRegressor follows the same struct layout as OLSRegressor: a fitted flag, nFeatures, intercept, and a coef vector. Two extra diagnostic fields, rank and condEst, are populated after every successful fit and are useful for inspecting the numerical health of the design.

struct QRRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    int rank = 0;       // p after a full-rank fit; may be less if pivoting were used
    double condEst = 0.0;  // max|R_kk| / min|R_kk|

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

The core of fit builds the augmented matrix, applies $p$ Householder reflections, checks the rank, then back-substitutes:

void fit(const Matrix& X, const vector<double>& y) {
    fitted = false;
    validateFitInput(X, y, "QRRegressor");
    int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;

    // Build X_b: intercept column first, then caller's features.
    Matrix Xb(n, vector<double>(p));
    for (int i = 0; i < n; ++i) {
        Xb[i][0] = 1.0;
        for (int j = 0; j < d; ++j) Xb[i][j+1] = X[i][j];
    }
    vector<double> rhs = y;   // accumulates Q^T y in-place

    const double eps = numeric_limits<double>::epsilon();
    for (int k = 0; k < p; ++k) {
        // Norm of sub-column Xb[k:n][k].
        double sigma2 = 0.0;
        for (int i = k; i < n; ++i) sigma2 += Xb[i][k] * Xb[i][k];
        if (!isfinite(sigma2))
            throw runtime_error(
                "QRRegressor::fit: column norm became non-finite; rescale the input");
        double sigma = sqrt(sigma2);

        if (sigma < eps * 64.0)
            throw runtime_error(
                "QRRegressor::fit: column " + to_string(k) +
                " is linearly dependent; design is rank-deficient. "
                "Use SVDRegressor for a min-norm solution.");

        // Reflector v: choose sign to avoid cancellation.
        double alpha = -copysign(sigma, Xb[k][k]);
        vector<double> v(n - k);
        for (int i = k; i < n; ++i) v[i-k] = Xb[i][k];
        v[0] -= alpha;  // v[0] = Xb[k][k] - alpha

        double vTv = 0.0;
        for (double vi : v) vTv += vi * vi;
        if (!isfinite(vTv))
            throw runtime_error(
                "QRRegressor::fit: reflector norm became non-finite; rescale the input");
        if (vTv < eps * eps) continue;

        // Apply H = I - 2vv^T/v^Tv to columns k..p-1 of Xb.
        for (int j = k; j < p; ++j) {
            double vTx = 0.0;
            for (int i = k; i < n; ++i) vTx += v[i-k] * Xb[i][j];
            double tau = 2.0 * vTx / vTv;
            for (int i = k; i < n; ++i) Xb[i][j] -= tau * v[i-k];
        }
        // Apply same reflector to accumulated Q^T y.
        double vTr = 0.0;
        for (int i = k; i < n; ++i) vTr += v[i-k] * rhs[i];
        double tau = 2.0 * vTr / vTv;
        for (int i = k; i < n; ++i) rhs[i] -= tau * v[i-k];
    }

    // Rank check via R diagonal.
    double maxAbs = 0.0;
    for (int k = 0; k < p; ++k) maxAbs = max(maxAbs, fabs(Xb[k][k]));
    double rankTol = maxAbs * (double)p * eps * 64.0;
    rank = 0;
    for (int k = 0; k < p; ++k) if (fabs(Xb[k][k]) > rankTol) ++rank;

    if (rank < p)
        throw runtime_error(
            "QRRegressor::fit: " + to_string(p - rank) +
            " R diagonal entr" + (p - rank == 1 ? "y" : "ies") +
            " below rank tolerance; design is numerically rank-deficient. "
            "Use SVDRegressor for a min-norm solution.");

    double minAbs = fabs(Xb[0][0]);
    for (int k = 0; k < p; ++k) minAbs = min(minAbs, fabs(Xb[k][k]));
    condEst = (minAbs > 0.0) ? maxAbs / minAbs
                              : numeric_limits<double>::infinity();

    // Back-substitute R * beta = (Q^T y)[0:p].
    vector<double> beta(p, 0.0);
    for (int k = p-1; k >= 0; --k) {
        double s = rhs[k];
        for (int j = k+1; j < p; ++j) s -= Xb[k][j] * beta[j];
        beta[k] = s / Xb[k][k];
        if (!isfinite(beta[k]))
            throw runtime_error(
                "QRRegressor::fit: back-substitution produced a non-finite coefficient");
    }

    intercept = beta[0];
    coef.assign(beta.begin()+1, beta.end());
    nFeatures = d;
    fitted = true;
}

One-Sided Jacobi SVD

The singular value decomposition $X_b = USV^\top$ factors the $n \times p$ augmented matrix into an orthogonal $U$ (columns are the left singular vectors), a diagonal $S$ with nonnegative entries (singular values in some order), and an orthogonal $V$ (columns are the right singular vectors). The least-squares solution, including a minimum-norm solution for rank-deficient problems, is the pseudoinverse:

where $S^+$ replaces each nonzero singular value with its reciprocal and leaves the zero (or near-zero) singular values at zero. This truncated form is the minimum-norm least-squares solution: among all $\beta$ that minimize $\lVert y - X_b\beta \rVert_2$, it selects the one with the smallest $\lVert \beta \rVert_2$. Singular values below a tolerance are omitted from the sum, and their corresponding directions in $\mathbb{R}^p$ (the null space of $X_b$) are not constrained by the data.

The Jacobi Iteration on Columns of $X_b$

One-sided Jacobi SVD never forms $X_b^\top X_b$. Instead it applies plane rotations directly to the columns of a working copy of $X_b$ and accumulates the same rotations in $V$, driving the columns of the working matrix toward orthogonality. Once converged, the column norms equal the singular values and the normalized columns are the left singular vectors $U$.

At each iteration step the algorithm selects a pair of columns $(i, j)$ and computes the inner products of the Gram matrix fragment: $G_{ii} = \lVert a_i \rVert^2$, $G_{jj} = \lVert a_j \rVert^2$, and $G_{ij} = a_i^\top a_j$ where $a_i, a_j$ are columns $i$ and $j$ of the current working matrix. If $|G_{ij}| / \sqrt{G_{ii} G_{jj}}$ is already below a tight tolerance, the pair is skipped as already orthogonal. Otherwise the Jacobi rotation angle is

and the rotation is applied to both the working matrix and $V$:

The algorithm sweeps over all $\binom{p}{2}$ column pairs, repeating until no pair triggers a rotation. Convergence is typically reached in a handful of sweeps for the small $p$ values common in regression problems. Each pair costs $O(n)$ to measure and rotate the working columns plus $O(p)$ to rotate $V$, so one full sweep costs $O(np^2 + p^3)$. In the tall-matrix regime $n \ge p$, this is $O(np^2)$, and $s$ sweeps cost $O(snp^2)$.

Rank, Tolerance, and Min-Norm

After convergence the singular values are extracted as $\sigma_k = \lVert a_k \rVert_2$ (the column norms of the working matrix), then sorted in descending order along with the corresponding columns of $V$. A rank tolerance is set as

where $\sigma_0$ is the largest singular value. Singular values below this threshold are treated as numerically zero, and the effective rank is the count of values above it. The condition estimate is $\sigma_0 / \sigma_{r-1}$ where $r$ is the effective rank. For a full-rank design this equals the true condition number of $X_b$; for a rank-deficient design it measures the conditioning within the identifiable subspace.

The regression coefficients are then

where $a_k^\top y / \sigma_k$ plays the role of $u_k^\top y$ (the normalized dot product) and the extra $1/\sigma_k$ factor comes from the pseudoinverse. The first element of $\beta$ is the intercept; the remaining $d$ elements are the feature coefficients.

One Sweep by Hand

The Jacobi rotation is easiest to trust after watching a single pair go orthogonal. A compact augmented matrix with deliberately correlated columns is $A = \begin{bmatrix}1&2\\1&3\\1&4\end{bmatrix}$: column $a_0 = [1,1,1]^\top$ is the intercept, and column $a_1 = [2,3,4]^\top$ is a feature whose nonzero mean makes it strongly aligned with the intercept. With $p = 2$ there is exactly one column pair $(0,1)$, so one rotation completes the sweep.

1. Gram fragment. $G_{00} = \lVert a_0 \rVert^2 = 3$, $G_{11} = \lVert a_1 \rVert^2 = 4+9+16 = 29$, and $G_{01} = a_0^\top a_1 = 2+3+4 = 9$. These three inner products come straight from the working columns; the full matrix $A^\top A$ is never assembled.
2. Normalized correlation. $|G_{01}| / \sqrt{G_{00} G_{11}} = 9 / \sqrt{3 \cdot 29} = 9/\sqrt{87} \approx 0.9649$. This is far above the skip threshold $4\,\epsilon_{\text{mach}}$, so the pair is rotated rather than left alone.
3. Stable rotation. $\zeta = (G_{11}-G_{00}) / (2 G_{01}) = 26/18 \approx 1.4444$. The tangent uses the sign-preserving, cancellation-free form $t = \operatorname{sign}(\zeta) / (|\zeta| + \sqrt{1+\zeta^2}) \approx 0.3124$, giving $c = 1/\sqrt{1+t^2} \approx 0.9545$ and $s = tc \approx 0.2982$.
4. Orthogonal columns. Rotating $a_0 \leftarrow c\,a_0 - s\,a_1$ and $a_1 \leftarrow s\,a_0 + c\,a_1$ yields new columns $[0.3582,\,0.0600,\,-0.2382]^\top$ and $[2.2072,\,3.1617,\,4.1162]^\top$. Their inner product is $\approx 1.2 \times 10^{-15}$, numerically zero, so no further rotation of this pair is needed.
5. Singular values and sorting. The new column norms are $0.4343$ and $5.6402$; these are the singular values. Sorting in descending order swaps them, so $\sigma_0 = 5.6402$ (the former column 1) and $\sigma_1 = 0.4343$, and the columns of the accumulated $V$ are reordered to match.
6. Rank tolerance and contributing directions. With $\text{tol} = \sigma_0 \cdot \max(n,p) \cdot \epsilon_{\text{mach}} \cdot 4 \approx 1.5 \times 10^{-14}$, both singular values exceed the threshold, so the effective rank is $2$ and the condition estimate is $\sigma_0/\sigma_1 \approx 12.99$. Both singular directions clear the tolerance, so both terms $\tfrac{a_k^\top y}{\sigma_k^2}\,v_k$ contribute to the pseudoinverse sum; nothing is truncated for this full-rank design.

The animation replays this exact sweep: Gram measurement, the stable rotation, the columns snapping orthogonal, and the descending sort that fixes the singular-value order feeding the pseudoinverse.

SVD Regressor Implementation

struct SVDRegressor {
    bool fitted = false;
    int nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    int rank = 0;
    double condEst = 0.0;
    double singularTol = 0.0;  // 0 = auto from data scale

    void fit(const Matrix& X, const vector<double>& y) {
        fitted = false;
        validateFitInput(X, y, "SVDRegressor");
        if (!isfinite(singularTol) || singularTol < 0.0)
            throw invalid_argument(
                "SVDRegressor::fit: singularTol must be finite and non-negative");
        int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;

        Matrix A(n, vector<double>(p));
        for (int i = 0; i < n; ++i) {
            A[i][0] = 1.0;
            for (int j = 0; j < d; ++j) A[i][j+1] = X[i][j];
        }
        Matrix V(p, vector<double>(p, 0.0));
        for (int k = 0; k < p; ++k) V[k][k] = 1.0;  // V = I_p

        const double eps = numeric_limits<double>::epsilon();
        const int maxSweeps = 100;
        bool changed = true;
        int sweep = 0;
        for (; sweep < maxSweeps && changed; ++sweep) {
            changed = false;
            for (int i = 0; i < p-1; ++i) {
                for (int j = i+1; j < p; ++j) {
                    double Gii = 0.0, Gjj = 0.0, Gij = 0.0;
                    for (int row = 0; row < n; ++row) {
                        Gii += A[row][i] * A[row][i];
                        Gjj += A[row][j] * A[row][j];
                        Gij += A[row][i] * A[row][j];
                    }
                    if (!isfinite(Gii) || !isfinite(Gjj) || !isfinite(Gij))
                        throw runtime_error(
                            "SVDRegressor::fit: Jacobi inner product became non-finite; "
                            "rescale the input");
                    double denom = sqrt(Gii) * sqrt(Gjj);
                    if (denom == 0.0 || fabs(Gij) / denom < eps * 4.0) continue;

                    double zeta = (Gjj - Gii) / (2.0 * Gij);
                    double azeta = fabs(zeta), t;
                    if (azeta > 1.0) {
                        double inv = 1.0 / azeta;
                        t = copysign(inv / (1.0 + hypot(1.0, inv)), zeta);
                    } else {
                        t = copysign(1.0 / (azeta + hypot(1.0, azeta)), zeta);
                    }
                    if (t == 0.0) continue;
                    double c = 1.0 / sqrt(1.0 + t*t);
                    double s = t * c;

                    for (int row = 0; row < n; ++row) {
                        double ai = A[row][i], aj = A[row][j];
                        A[row][i] = c*ai - s*aj;
                        A[row][j] = s*ai + c*aj;
                    }
                    for (int row = 0; row < p; ++row) {
                        double vi = V[row][i], vj = V[row][j];
                        V[row][i] = c*vi - s*vj;
                        V[row][j] = s*vi + c*vj;
                    }
                    changed = true;
                }
            }
        }
        if (changed)
            throw runtime_error(
                "SVDRegressor::fit: Jacobi iteration did not converge in " +
                to_string(maxSweeps) + " sweeps");

        // Extract and sort singular values.
        vector<double> S(p);
        for (int k = 0; k < p; ++k) {
            double s2 = 0.0;
            for (int row = 0; row < n; ++row) s2 += A[row][k] * A[row][k];
            if (!isfinite(s2))
                throw runtime_error(
                    "SVDRegressor::fit: singular-value norm became non-finite");
            S[k] = sqrt(s2);
        }
        vector<int> order(p);
        iota(order.begin(), order.end(), 0);
        sort(order.begin(), order.end(), [&](int a, int b){ return S[a] > S[b]; });

        vector<double> Ss(p);
        Matrix As(n, vector<double>(p));
        Matrix Vs(p, vector<double>(p));
        for (int k = 0; k < p; ++k) {
            Ss[k] = S[order[k]];
            for (int row = 0; row < n; ++row) As[row][k] = A[row][order[k]];
            for (int row = 0; row < p; ++row) Vs[row][k] = V[row][order[k]];
        }

        double tol = singularTol > 0.0
                     ? singularTol
                     : Ss[0] * (double)max(n, p) * eps * 4.0;
        rank = 0;
        for (int k = 0; k < p; ++k) if (Ss[k] > tol) ++rank;

        if (rank == 0)
            throw runtime_error(
                "SVDRegressor::fit: all singular values are below tolerance; "
                "check input for near-zero or identical features");

        condEst = Ss[0] / Ss[rank-1];

        // beta = sum_{k < rank} (A[:,k] . y / S_k^2) * V[:,k]
        vector<double> beta(p, 0.0);
        for (int k = 0; k < rank; ++k) {
            double coeff = 0.0;
            for (int row = 0; row < n; ++row) coeff += As[row][k] * y[row];
            coeff /= (Ss[k] * Ss[k]);
            for (int row = 0; row < p; ++row) beta[row] += coeff * Vs[row][k];
        }
        for (int k = 0; k < p; ++k)
            if (!isfinite(beta[k]))
                throw runtime_error(
                    "SVDRegressor::fit: pseudoinverse produced a non-finite coefficient");

        intercept = beta[0];
        coef.assign(beta.begin()+1, beta.end());
        nFeatures = d;
        fitted = true;
    }
    // predict overloads identical to QRRegressor
};

Validation and Errors

Both solvers share a validation helper that checks input shapes and finiteness before any computation begins, matching the convention from the OLS post:

• Row count mismatch between X and y throws std::invalid_argument.
• An empty training set, zero feature columns, ragged rows, or non-finite entries all throw std::invalid_argument with the row index named where relevant.
• Fewer than $p = d{+}1$ rows throws std::invalid_argument before any reflector or rotation is computed.
• QRRegressor throws std::runtime_error if any Householder sub-column norm is numerically zero (immediately rank-deficient) or if the post-factorization rank check detects near-linear dependence.
• SVDRegressor rejects a negative or non-finite custom tolerance and throws std::runtime_error if Jacobi iteration misses its sweep cap or all singular values fall below tolerance.
• Calling predict before fit throws std::logic_error; a feature-count mismatch or non-finite feature throws std::invalid_argument.

Three-Scenario Comparison

The complete program below runs both solvers against three designs on the same synthetic true model $y = 3 + 2x_1 - 1.5x_2 + \varepsilon$ (same seed and split as the OLS post). Scenario A is well-conditioned independent features. Scenario B has $x_2 \approx x_1$, introducing near-collinearity. Scenario C has $x_2 = 2x_1$ exactly, making the design rank-deficient. The output shows that QR and SVD agree with OLS on Scenario A, both flag the condition inflation in Scenario B, QR throws on Scenario C while SVD computes the minimum-norm solution.

#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>>;

// ============================================================
// Input validation helper shared by both solvers.
// ============================================================
static void validateFitInput(const Matrix& X, const vector<double>& y,
                              const string& who) {
    if (X.size() != y.size())
        throw invalid_argument(who + "::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(who + "::fit: empty training set");
    int d = (int)X[0].size();
    if (d == 0) throw invalid_argument(who + "::fit: zero feature columns");
    for (int i = 0; i < n; ++i) {
        if ((int)X[i].size() != d)
            throw invalid_argument(who + "::fit: row " + to_string(i) +
                                   " has " + to_string(X[i].size()) +
                                   " columns, expected " + to_string(d));
        for (int j = 0; j < d; ++j)
            if (!isfinite(X[i][j]))
                throw invalid_argument(who + "::fit: non-finite feature in row " + to_string(i));
        if (!isfinite(y[i]))
            throw invalid_argument(who + "::fit: non-finite target at row " + to_string(i));
    }
    if (n < d + 1)
        throw invalid_argument(who + "::fit: need at least " + to_string(d + 1) +
                               " rows for " + to_string(d) + " features plus intercept, got " +
                               to_string(n));
}

// ============================================================
// Metric helpers.
// ============================================================
static 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));
}

static 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();
}

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

static 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;
}

// ============================================================
// QRRegressor: Householder triangularization on X_b directly.
// ============================================================
struct QRRegressor {
    bool   fitted   = false;
    int    nFeatures = 0;
    double intercept = 0.0;
    vector<double> coef;
    int    rank    = 0;
    double condEst = 0.0;

    void fit(const Matrix& X, const vector<double>& y) {
        fitted = false;
        validateFitInput(X, y, "QRRegressor");
        int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;

        // Build augmented matrix X_b (intercept column first).
        Matrix Xb(n, vector<double>(p));
        for (int i = 0; i < n; ++i) {
            Xb[i][0] = 1.0;
            for (int j = 0; j < d; ++j) Xb[i][j+1] = X[i][j];
        }
        vector<double> rhs = y;  // accumulates Q^T y in-place

        const double eps = numeric_limits<double>::epsilon();

        for (int k = 0; k < p; ++k) {
            // Compute norm of sub-column Xb[k:n][k].
            double sigma2 = 0.0;
            for (int i = k; i < n; ++i) sigma2 += Xb[i][k] * Xb[i][k];
            if (!isfinite(sigma2))
                throw runtime_error(
                    "QRRegressor::fit: column norm became non-finite; rescale the input");
            double sigma = sqrt(sigma2);

            if (sigma < eps * 64.0)
                throw runtime_error(
                    "QRRegressor::fit: column " + to_string(k) +
                    " is linearly dependent; design is rank-deficient. "
                    "Use SVDRegressor for a min-norm solution.");

            // Householder reflector: choose sign to avoid cancellation.
            double alpha = -copysign(sigma, Xb[k][k]);
            vector<double> v(n - k);
            for (int i = k; i < n; ++i) v[i-k] = Xb[i][k];
            v[0] -= alpha;

            double vTv = 0.0;
            for (double vi : v) vTv += vi * vi;
            if (!isfinite(vTv))
                throw runtime_error(
                    "QRRegressor::fit: reflector norm became non-finite; rescale the input");
            if (vTv < eps * eps) continue;

            // Apply H = I - 2vv^T/v^Tv to columns k..p-1 of Xb.
            for (int j = k; j < p; ++j) {
                double vTx = 0.0;
                for (int i = k; i < n; ++i) vTx += v[i-k] * Xb[i][j];
                double tau = 2.0 * vTx / vTv;
                for (int i = k; i < n; ++i) Xb[i][j] -= tau * v[i-k];
            }
            // Apply same reflector to rhs (accumulates Q^T y).
            double vTr = 0.0;
            for (int i = k; i < n; ++i) vTr += v[i-k] * rhs[i];
            double tau = 2.0 * vTr / vTv;
            for (int i = k; i < n; ++i) rhs[i] -= tau * v[i-k];
        }

        // Rank check via R diagonal magnitudes.
        double maxAbs = 0.0;
        for (int k = 0; k < p; ++k) maxAbs = max(maxAbs, fabs(Xb[k][k]));
        double rankTol = maxAbs * (double)p * eps * 64.0;
        rank = 0;
        for (int k = 0; k < p; ++k) if (fabs(Xb[k][k]) > rankTol) ++rank;

        if (rank < p)
            throw runtime_error(
                "QRRegressor::fit: " + to_string(p - rank) +
                " R diagonal entr" + (p - rank == 1 ? "y" : "ies") +
                " below rank tolerance; design is numerically rank-deficient. "
                "Use SVDRegressor for a min-norm solution.");

        double minAbs = fabs(Xb[0][0]);
        for (int k = 1; k < p; ++k) minAbs = min(minAbs, fabs(Xb[k][k]));
        condEst = (minAbs > 0.0) ? maxAbs / minAbs
                                  : numeric_limits<double>::infinity();

        // Back-substitute R * beta = (Q^T y)[0:p].
        vector<double> beta(p, 0.0);
        for (int k = p-1; k >= 0; --k) {
            double s = rhs[k];
            for (int j = k+1; j < p; ++j) s -= Xb[k][j] * beta[j];
            beta[k] = s / Xb[k][k];
            if (!isfinite(beta[k]))
                throw runtime_error(
                    "QRRegressor::fit: back-substitution produced a non-finite coefficient");
        }

        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("QRRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("QRRegressor::predict: expected " +
                                   to_string(nFeatures) + " features, got " +
                                   to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j]))
                throw invalid_argument(
                    "QRRegressor::predict: non-finite feature at index " +
                    to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out))
            throw runtime_error("QRRegressor::predict: prediction became non-finite");
        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;
    }
};

// ============================================================
// SVDRegressor: one-sided Jacobi SVD directly on X_b.
// ============================================================
struct SVDRegressor {
    bool   fitted      = false;
    int    nFeatures   = 0;
    double intercept   = 0.0;
    vector<double> coef;
    int    rank        = 0;
    double condEst     = 0.0;
    double singularTol = 0.0;  // 0 = auto from data scale

    void fit(const Matrix& X, const vector<double>& y) {
        fitted = false;
        validateFitInput(X, y, "SVDRegressor");
        if (!isfinite(singularTol) || singularTol < 0.0)
            throw invalid_argument(
                "SVDRegressor::fit: singularTol must be finite and non-negative");
        int n = (int)X.size(), d = (int)X[0].size(), p = d + 1;

        // Build working copy of X_b.
        Matrix A(n, vector<double>(p));
        for (int i = 0; i < n; ++i) {
            A[i][0] = 1.0;
            for (int j = 0; j < d; ++j) A[i][j+1] = X[i][j];
        }
        // V starts as identity.
        Matrix V(p, vector<double>(p, 0.0));
        for (int k = 0; k < p; ++k) V[k][k] = 1.0;

        const double eps = numeric_limits<double>::epsilon();
        const int maxSweeps = 100;
        bool changed = true;
        int sweep = 0;
        for (; sweep < maxSweeps && changed; ++sweep) {
            changed = false;
            for (int i = 0; i < p-1; ++i) {
                for (int j = i+1; j < p; ++j) {
                    double Gii = 0.0, Gjj = 0.0, Gij = 0.0;
                    for (int row = 0; row < n; ++row) {
                        Gii += A[row][i] * A[row][i];
                        Gjj += A[row][j] * A[row][j];
                        Gij += A[row][i] * A[row][j];
                    }
                    if (!isfinite(Gii) || !isfinite(Gjj) || !isfinite(Gij))
                        throw runtime_error(
                            "SVDRegressor::fit: Jacobi inner product became non-finite; "
                            "rescale the input");
                    double denom = sqrt(Gii) * sqrt(Gjj);
                    if (denom == 0.0 || fabs(Gij) / denom < eps * 4.0) continue;

                    double zeta = (Gjj - Gii) / (2.0 * Gij);
                    double azeta = fabs(zeta), t;
                    if (azeta > 1.0) {
                        double inv = 1.0 / azeta;
                        t = copysign(inv / (1.0 + hypot(1.0, inv)), zeta);
                    } else {
                        t = copysign(1.0 / (azeta + hypot(1.0, azeta)), zeta);
                    }
                    if (t == 0.0) continue;
                    double c = 1.0 / sqrt(1.0 + t*t);
                    double s = t * c;

                    for (int row = 0; row < n; ++row) {
                        double ai = A[row][i], aj = A[row][j];
                        A[row][i] = c*ai - s*aj;
                        A[row][j] = s*ai + c*aj;
                    }
                    for (int row = 0; row < p; ++row) {
                        double vi = V[row][i], vj = V[row][j];
                        V[row][i] = c*vi - s*vj;
                        V[row][j] = s*vi + c*vj;
                    }
                    changed = true;
                }
            }
        }
        if (changed)
            throw runtime_error(
                "SVDRegressor::fit: Jacobi iteration did not converge in " +
                to_string(maxSweeps) + " sweeps");

        // Extract singular values (column norms of A after convergence).
        vector<double> S(p);
        for (int k = 0; k < p; ++k) {
            double s2 = 0.0;
            for (int row = 0; row < n; ++row) s2 += A[row][k] * A[row][k];
            if (!isfinite(s2))
                throw runtime_error(
                    "SVDRegressor::fit: singular-value norm became non-finite");
            S[k] = sqrt(s2);
        }

        // Sort singular triplets in descending order.
        vector<int> order(p);
        iota(order.begin(), order.end(), 0);
        sort(order.begin(), order.end(), [&](int a, int b){ return S[a] > S[b]; });

        vector<double> Ss(p);
        Matrix As(n, vector<double>(p));
        Matrix Vs(p, vector<double>(p));
        for (int k = 0; k < p; ++k) {
            Ss[k] = S[order[k]];
            for (int row = 0; row < n; ++row) As[row][k] = A[row][order[k]];
            for (int row = 0; row < p; ++row) Vs[row][k] = V[row][order[k]];
        }

        double tol = singularTol > 0.0
                     ? singularTol
                     : Ss[0] * (double)max(n, p) * eps * 4.0;
        rank = 0;
        for (int k = 0; k < p; ++k) if (Ss[k] > tol) ++rank;

        if (rank == 0)
            throw runtime_error(
                "SVDRegressor::fit: all singular values are below tolerance; "
                "check input for near-zero or identical features");

        condEst = Ss[0] / Ss[rank-1];

        // beta = sum_{k < rank} (As[:,k] . y / Ss[k]^2) * Vs[:,k]
        vector<double> beta(p, 0.0);
        for (int k = 0; k < rank; ++k) {
            double coeff = 0.0;
            for (int row = 0; row < n; ++row) coeff += As[row][k] * y[row];
            coeff /= (Ss[k] * Ss[k]);
            for (int row = 0; row < p; ++row) beta[row] += coeff * Vs[row][k];
        }
        for (int k = 0; k < p; ++k)
            if (!isfinite(beta[k]))
                throw runtime_error(
                    "SVDRegressor::fit: pseudoinverse produced a non-finite coefficient");

        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("SVDRegressor::predict: model is not fitted");
        if ((int)x.size() != nFeatures)
            throw invalid_argument("SVDRegressor::predict: expected " +
                                   to_string(nFeatures) + " features, got " +
                                   to_string(x.size()));
        double out = intercept;
        for (int j = 0; j < nFeatures; ++j) {
            if (!isfinite(x[j]))
                throw invalid_argument(
                    "SVDRegressor::predict: non-finite feature at index " +
                    to_string(j));
            out += coef[j] * x[j];
        }
        if (!isfinite(out))
            throw runtime_error("SVDRegressor::predict: prediction became non-finite");
        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;
    }
};

// ============================================================
// main: three scenarios
// ============================================================
int main() {
    // === Scenario A: well-conditioned (2 independent features) ===
    printf("=== Scenario A: well-conditioned (2 independent features) ===\n");
    {
        mt19937 rng(7);
        uniform_real_distribution<double> ux(-3.0, 3.0);
        normal_distribution<double> noise(0.0, 0.6);
        const int N = 240;
        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(3.0 + 2.0*x1 - 1.5*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]); }
        }
        // QR
        QRRegressor qr; qr.fit(Xtr, ytr);
        auto ptQR = qr.predict(Xte);
        printf("[QR]  intercept=%.4f  coef=[%.4f, %.4f]  cond=%.1f  rank=%d\n",
               qr.intercept, qr.coef[0], qr.coef[1], qr.condEst, qr.rank);
        printf("[QR]  test RMSE=%.4f  test R^2=%.4f\n", rmse(yte, ptQR), r2(yte, ptQR));
        // SVD
        SVDRegressor sv; sv.fit(Xtr, ytr);
        auto ptSV = sv.predict(Xte);
        printf("[SVD] intercept=%.4f  coef=[%.4f, %.4f]  cond=%.1f  rank=%d\n",
               sv.intercept, sv.coef[0], sv.coef[1], sv.condEst, sv.rank);
        printf("[SVD] test RMSE=%.4f  test R^2=%.4f\n", rmse(yte, ptSV), r2(yte, ptSV));
    }
    // === Scenario B: nearly collinear features (x2 ~ x1) ===
    printf("\n=== Scenario B: nearly collinear features (x2 ~ x1) ===\n");
    {
        mt19937 rng(42);
        normal_distribution<double> base(0.0, 1.0);
        normal_distribution<double> perturb(0.0, 0.02);
        normal_distribution<double> yNoise(0.0, 0.3);
        const int N = 200;
        Matrix X; vector<double> y;
        for (int i = 0; i < N; ++i) {
            double x1 = base(rng), x2 = x1 + perturb(rng);
            X.push_back({x1, x2});
            y.push_back(1.0 + 2.0*x1 - 0.5*x2 + yNoise(rng));
        }
        QRRegressor qr; qr.fit(X, y);
        auto pQR = qr.predict(X);
        printf("[QR]  intercept=%.4f  coef=[%.4f, %.4f]  cond=%.1f  rank=%d\n",
               qr.intercept, qr.coef[0], qr.coef[1], qr.condEst, qr.rank);
        printf("[QR]  train RMSE=%.4f  R^2=%.4f\n", rmse(y, pQR), r2(y, pQR));
        SVDRegressor sv; sv.fit(X, y);
        auto pSV = sv.predict(X);
        printf("[SVD] intercept=%.4f  coef=[%.4f, %.4f]  cond=%.1f  rank=%d\n",
               sv.intercept, sv.coef[0], sv.coef[1], sv.condEst, sv.rank);
        printf("[SVD] train RMSE=%.4f  R^2=%.4f\n", rmse(y, pSV), r2(y, pSV));
    }
    // === Scenario C: exact rank deficiency (x2 = 2*x1) ===
    printf("\n=== Scenario C: exact rank deficiency (x2 = 2*x1) ===\n");
    {
        mt19937 rng(99);
        normal_distribution<double> base(0.0, 1.0);
        normal_distribution<double> yNoise(0.0, 0.3);
        const int N = 100;
        Matrix X; vector<double> y;
        for (int i = 0; i < N; ++i) {
            double x1 = base(rng);
            X.push_back({x1, 2.0*x1});        // exactly proportional
            y.push_back(2.0 + 1.5*x1 - 1.0*(2.0*x1) + yNoise(rng));
        }
        try {
            QRRegressor qr; qr.fit(X, y);
        } catch (const runtime_error& e) {
            printf("[QR]  caught: %s\n", e.what());
        }
        SVDRegressor sv; sv.fit(X, y);
        auto pSV = sv.predict(X);
        printf("[SVD] intercept=%.4f  coef=[%.4f, %.4f]  rank=%d  cond=%.2e\n",
               sv.intercept, sv.coef[0], sv.coef[1], sv.rank, sv.condEst);
        printf("[SVD] coef[0]+2*coef[1]=%.4f (true -0.5)  RMSE=%.4f  R^2=%.4f\n",
               sv.coef[0] + 2.0*sv.coef[1], rmse(y, pSV), r2(y, pSV));
    }
    return 0;
}

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

Scenario A shows QR and SVD recovering the same coefficients as OLS: intercept near $3.0$, $\beta_1 \approx 2.0$, $\beta_2 \approx -1.5$, both with the same test metrics. The condition estimate of $1.8$ confirms this design is very well-conditioned. Scenario B recovers a good fit ($R^2 = 0.96$) but the individual coefficients are far from the true values ($2.0$ and $-0.5$) because with $x_2 \approx x_1$ the data cannot distinguish their separate effects. The QR condition estimate of $47.1$ and the SVD estimate of $93.7$ reflect different measures: QR reports the ratio of the largest to smallest $|R_{kk}|$ diagonal entry, which is a lower bound; SVD reports $\sigma_{\max}/\sigma_{\min}$, the exact condition of $X_b$. Scenario C has $x_2 = 2x_1$, so the design is exactly rank-2 out of a possible rank-3. QR throws at the point where the zero sub-column is detected. SVD truncates at the near-zero singular value and recovers $\beta_1 + 2\beta_2 \approx -0.49$ (true $-0.5$), which is the only identifiable combination; the individual values $[-0.097, -0.195]$ are the minimum-norm split of that identifiable direction.

The Minimum-Norm Split by Hand

Scenario C rewards a slower look, because it shows exactly what "minimum norm" buys when a design is rank-deficient. The feature matrix has $x_2 = 2x_1$, so the augmented columns are $[\mathbf{1},\, x_1,\, 2x_1]$. Any prediction depends on the features only through the combination $\beta_1 x_1 + \beta_2 (2x_1) = (\beta_1 + 2\beta_2)\,x_1$, so the single number $g = \beta_1 + 2\beta_2$ is the only feature-slope quantity the data can identify. The output reports $g \approx -0.4869$, close to the true $-0.5$.

That leaves one degree of freedom completely unconstrained by the data: the null space. A vector $z$ lies in the null space of $X_b$ when $X_b z = 0$ for every row, which forces $z_0 = 0$ and $z_1 + 2z_2 = 0$, i.e. $z \propto (0,\,-2,\,1)$. Moving $\beta$ along this direction changes the individual coefficients while leaving every prediction, and therefore the loss, untouched. The pseudoinverse resolves the ambiguity by selecting the solution orthogonal to the null space, which is the minimum-norm choice:

Substituting the null-space constraint into the identifiable combination pins down both coefficients from the single measured number:

Those are exactly the printed coefficients $[-0.0974, -0.1948]$, and the ratio $\beta_2/\beta_1 = 2$ confirms the solution sits on the row space rather than drifting along the null space. A different but equally valid least-squares answer, for example $\beta_1 = -0.4869,\ \beta_2 = 0$, would predict identically yet carry a larger Euclidean norm; the pseudoinverse rejects it in favor of the shortest vector consistent with the data. The intercept stays separately identifiable at $\beta_0 \approx 2.032$ because the ones column is independent of $x_1$, so only the two feature coefficients participate in the split.

When to Use Each Solver

Complexity and Memory

Both solvers are $O(np^2)$ in the regime $n \gg p$, which matches the normal-equation path once the $O(nd^2)$ assembly step is included. The extra $O(np)$ allocation for $X_b$ is the main additional cost over the OLS approach, which avoided materializing $X_b$ at all.

Practical Diagnostics

Both solvers expose enough state to sanity-check a fit without leaving the two implementations. Every check below uses only fields already computed by QRRegressor or SVDRegressor and the shared predict path.

• Residual orthogonality. A correct least-squares solution leaves a residual $r = y - X_b\hat\beta$ that is orthogonal to the column space of $X_b$, so $X_b^\top r \approx 0$. Concretely, form $r_i = y_i - \texttt{predict}(x_i)$ over the training rows and confirm that $\sum_i r_i \approx 0$ (orthogonality to the intercept column) and $\sum_i r_i x_{ij} \approx 0$ for each feature $j$ (orthogonality to feature column $j$). On Scenario A both sums land within a few multiples of $\epsilon_{\text{mach}}$ times the data scale; a residual that correlates with a feature column signals a solve that stopped short of the true minimizer.
• Reconstruction and column orthogonality. The one-sided Jacobi routine reaches convergence only when the working columns are mutually orthogonal, so the natural SVD check is $a_i^\top a_j \approx 0$ for $i \ne j$ after the sweeps, which the hand walkthrough verified at $\approx 10^{-15}$. Because the routine keeps the accumulated $V$ and the rotated columns $A_{\text{final}}$, the identity $A_{\text{final}} = X_b V$ reconstructs the input design from the singular triplets. For QR the analogous reconstruction is $R = Q^\top X_b$: the reflected feature column and right-hand side computed by hand are precisely the rows of that product, and a nonzero entry below the diagonal of the stored $R$ would flag a reflection that failed to zero its sub-column.
• QR rough condition versus SVD singular estimate. The two condEst fields are not measuring the same quantity. QR reports $\max_k|R_{kk}| / \min_k|R_{kk}|$, a cheap lower bound that only sees the diagonal of $R$; Scenario B prints $47.1$. SVD reports $\sigma_0/\sigma_{r-1}$, the exact ratio of extreme singular values of $X_b$; Scenario B prints $93.7$. When the two disagree by a wide margin the QR diagonal is understating the true conditioning, which is the expected behavior of an unpivoted factorization and a cue to trust the SVD figure.
• Tolerance sensitivity. The SVD rank is a thresholded count, so it is worth probing how stable that count is. Setting singularTol explicitly and sweeping it across a few decades reveals whether any singular value sits near the auto threshold $\sigma_0 \cdot \max(n,p)\cdot\epsilon_{\text{mach}}\cdot 4$. A rank that jumps as the tolerance crosses a value that is orders of magnitude below $\sigma_0$ indicates a genuinely well-separated spectrum; a rank that flips under a tiny tolerance change means the design is borderline rank-deficient and the reported coefficients in the smallest retained direction should be treated with caution.

Limitations and Edge Cases

• QR without pivoting and near-rank-deficiency. For near-zero but not exactly zero sub-columns, the Householder sigma check may not fire, and the rank check on $|R_{kk}|$ may also not fire if the small diagonal entry is above the tolerance. In this case the back-substitution will divide by a very small number and produce a coefficient vector with large, nearly-cancelling entries. Column-pivoting QR (not implemented here) reorders columns to bring the most informative directions first, guaranteeing that the largest available pivot is always used; this is the standard fix for near-rank-deficient designs in a QR context.
• SVD sweep count. For $p$ up to about $20$ the one-sided Jacobi algorithm typically converges in $5{-}10$ sweeps. For larger $p$ the convergence can slow, though it is always guaranteed for exact arithmetic. The cap of 100 sweeps is generous for any regression problem with double-precision data.
• Sign ambiguity in singular vectors. The sign of each column of $U$ and $V$ is arbitrary (flipping both leaves the decomposition valid). This does not affect the coefficients $\beta$ because the sign cancels in the pseudoinverse product, but it means the singular vector columns themselves should not be compared across runs without sign normalization.
• Min-norm vs other regularized solutions. The SVD minimum-norm solution is not the same as a ridge-regularized solution, which also produces a unique answer for rank-deficient designs but shrinks all coefficients toward zero rather than minimizing their Euclidean norm. The regularized regression post covers ridge and how it modifies the singular-value spectrum.
• Duplicate or constant feature columns. A column of all ones in the caller's $X$ duplicates the implicit intercept column; a constant feature column becomes linearly dependent on the intercept. Both cause rank deficiency. QR throws; SVD returns a min-norm fit but the resulting individual coefficients are not separately interpretable.

Next Step

Both QR and SVD are direct solvers: they process the entire training set in one factorization and produce an exact (or minimum-norm exact) solution. When $n$ or $d$ is large enough that forming and factoring an $n \times p$ matrix is expensive, repeated gradient updates provide a controllable alternative. The next post replaces the direct solve with an iterative approach. Its teaching batch fit retains a standardized training copy, while its frozen-scaler online path updates an existing model one sample at a time.