Why do OLS libraries fit models using the MP Pseudoinverse of the design matrix?

Question

For the linear model $y = X\beta$ for design matrix $X$, it's well known that the optimal solution is $\hat{\beta} = (X'X)^{-1}X'y$.

Some statistical libraries (such as Python's statsmodels) estimate parameters by first computing the pseudoinverse of the design matrix: $X^{+} = \text{pseudoinverse}(X) = (X'X)^{-1}X'$ and then computing the estimate as $\hat{\beta}=X^{+}y$.

I'm operating in a restricted memory environment with a "tall and skinny" design matrix with shapes up to (1e8,3). Computing the pseudoinverse requires allocating an additional matrix X^{+} of shape (3,1e8) before computing the final parameters and it is at this step that I often run out of memory.

We can contrast that with computing the parameters in 3 steps:

1. $A = (X'X)^{-1}$. This allocates a 3x3 matrix then inverts it.
2. $B = X'y$. This allocates a 3x1 vector.
3. $\hat{\beta} = AB$.

My question is: are there any downsides (such as numerical issues) to computing the OLS estimates using the methodology above? I can guarantee the absence of perfect multicollinearity (from the structure of the design matrix) thus guaranteeing that the inverse of the gram matrix ($A$) exists.

Answer 1

The libraries are making an assumption about what will make the most sense most of the time. The decision of the library is based on a memory-speed-numerical stability tradeoff that will make the most sense based on the shape and contents of $X$ the library expects to see. You have a very atypically shaped matrix $X$, so the library's assumption doesn't make sense for you. What you're proposing is perfectly valid and makes a lot of sense for your $X$ (especially if the columns aren't close to colinear).