I've started solving this problem.
Solve least squares problem X*theta_hat = y using Cholesky Decomposition.

|  Method   |   Operations    | Factor * np^2 |
|-----------|-----------------|---------------|
| Cholesky  |   1/3 * np^2    |      1/3      |
|    QR     |   p^3/3 + np^2  |   1 - p/3n    |
|    SVD    |   p^3   + np^2  |    1 - p/n    |

NOTE: HyperLearn's implementation of Cholesky Solve uses L2 Regularization to enforce stability.
Cholesky is known to fail on ill-conditioned problems, so adding L2 penalties helpes it.

Note, the algorithm in this implementation is as follows:

    alpha = dtype(X).decimal    [1e-6 is float32]
    while failure {
        solve cholesky ( XTX + alpha*identity )
        alpha *= 10
    }

If MSE (Mean Squared Error) is abnormally high, it might be better to solve using stabler but
slower methods like qr_solve, svd_solve or lstsq.

https://www.quora.com/Is-it-better-to-do-QR-Cholesky-or-SVD-for-solving-least-squares-estimate-and-why