[east2west]

Some clarifications. Explicitly computing matrix inversion is never a good idea, unless you require an inverted matrix. For solving the normal equation, indeed for solving any linear equation, use appropriate specific algorithms (variants of Gaussian elimination, QR, SVD, ...). For regression, computing the inverse of covariance matrix is actually faster than QR methods ( plural because there are more than one) at the cost of numerical instability. SVD is the more stable but also more expensive method. These are the standard methods you will find in any linear algebra or statistics textbooks. For regular use they are just fine. It is when you have thousands of variables and tens of thousands of measurements when they fall short. Iterative optimization algorithms can be used, or random matrix algorithms. There is no free lunch, because numerical instability is inherent in these algorithms, but you do have some, more relaxed, guarantee of errors or probability of errors.

The reason I said iterative optimization algorithms can be used for regression because linear regression is a quadratic optimization problem. We have analytic solutions because of special linear algebra properties. It can be said that we have a good handle on quadratic optimization problems, hence the multitude of algorithm choices. General optimization algorithms necessarily carry cost in numerical instability, but it may be a cost worth paying, especially for large data sets.

Finally, I want to echo your sentiment that DO NOT invert. Solve linear equations should be by tailored algorithms. In theory they yield identical results, in practice specific algorithms are much better.