Which is more numerically stable for OLS: pinv vs QR

If I am doing standard OLS and want to calculate beta values (OLS estimators), which of the following is the more numerically stable method? And why?

Assuming that the columns of $$X$$ are already mean-centered and normalised, to solve $$Y = X\beta + \epsilon$$ do:

1) $$\hat{\beta}_{pinv}=(X'X)^+X'Y$$

2) $$\hat{\beta}_{QR} = \text{solve}(R,Q'Y)$$

Where $$^+$$ represents the moore-penrose inverse, $$Q$$ and $$R$$ come from the QR decomposition of $$X$$ and solve is a function like the solve functions in python or r.

I would have thought (2) was better as $$(X'X)^+$$ seems to have a higher condition number than $$R$$, but in practice (in python at least) I am finding that the beta values derived from (1) minimize the sum of squared residuals better.

• When you write "minimize the sum of residuals" in your last line, do you mean "minimize the sum of squared residuals"? Commented Feb 27, 2019 at 19:27
• Apologies, yes I do. I will change this now. Commented Feb 27, 2019 at 20:32
• I cannot see why one method would "minimize the sum of squared residuals better"... In any case, QR is faster. QR-decomposition should be of order $\frac{4}{3} n^3$ while the SVD (the work-horse behind the MP-pseudo-inverse) takes $\frac{8}{3}n^3$. (For the record, Cholesky is the fastest with $\frac{1}{3}n^3$) You can quickly check RcppEigen's vignette for some timing informaton (check Table 2). Commented Feb 27, 2019 at 22:55
• When you say "computationally efficient" -- per your title and first paragraph -- you seem to be asking about calculation speed (perhaps measured in flops, or by multiplications and additions -- or maybe by some other measure), but the last paragraph of your question talks about condition number which seems to imply you're interested in something other than speed (in many situations one can get more speed by choosing a less stable algorithm). Can you be clearer about what you want to know? Commented Feb 27, 2019 at 23:11
• I am very sorry yes I was interested in the speed but my main concern was actually numerical stability - sorry to have caused such confusion. I have addressed this now. Thank you for the indepth analysis of the time in flops as well though - this is still greatly valued information and will go a long way in helping me with what I am doing. Commented Feb 28, 2019 at 1:53

Using the Moore-Penrose pseudo-inverse $$X^{\dagger}$$ of an matrix $$X$$ is more stable in the sense that can directly account for rank-deficient design matrices $$X$$. $$X^{\dagger}$$ allows us to naturally employ the identities: $$X^{\dagger} X X^{\dagger} = X$$ and $$X X^{\dagger} X= X^{\dagger}$$; the matrix $$X^{\dagger}$$ can be used as "surrogate" the true inverse of the matrix $$X$$, even if the inverse matrix $$X^{-1}$$ does not exist. In addition, the "usual" way of computing $$X^{\dagger}$$ by employing the Singular Value Decomposition of matrix $$X$$, where $$X = USV^T$$, is straight-forward methodologically and computationally well-studied. We simply take the reciprocal of the non-zero singular values in the diagonal matrix $$S$$, and we are good to go. Moore-Penrose pseudo-inverses are common in many proofs because they "just exist" and greatly simplify many derivations.
That said, in most cases it is not good practice to use the Moore-Penrose Pseudo-inverse unless we have a very good reason (e.g. our procedure consistently employs small and potentially rank-degenerate covariance matrices). The reasons are that: 1. It can hide true underlying problems with our data (e.g. duplication of variables) and 2. it is unnecessarily expensive (we have better alternatives). Finally, note that the Moore-Penrose pseudo-inverse of a full rank $$X$$ can be directed computed through the QR factorization of $$X$$, $$X = QR$$, as: $$X^{\dagger} = [R^{-1}_{1} 0] Q^T$$ where $$R_1$$ is an upper triangular matrix, coming from the "thin/reduced/skinny" QR factorization of $$X$$. So we do not really gain much if $$X$$ is full rank anyway. (Gentle's Matrix Algebra: Theory, Computations and Applications in Statistics provides a wealth of information the matter if one wishes to explore this further - Sect. 3.6 on Generalised Inverses should be a relevant starting point.)
To elaborate my first point a bit: It is far more natural to use a penalised regression procedure like Ridge or LASSO if we have issues with collinearity or simply have a $$p\gg n$$ (i.e. more predictors than data-points) than hide the problem using $$X^\dagger$$. The condition of a system of equations solved through the employment of Moore-Penrose pseudo-inverses might still be prohibitorily bad, resulting to unstable solutions and/or misleading inference. Thus if numerical stability is an issue, I would suggest using regularisation directly instead of Moore-Penrose pseudo-inverses. Note that in terms of speed, computing $$X^{\dagger}$$ is also problematic; potentially iterative methods based on gradient descent methods or alternating least squares are far faster for large systems (e.g. in Recommender Systems literature, see Paterek (2008) Improving regularized singular value decomposition for collaborative filtering for something very concise).