Add an additional OLS implementation based on Cholesky Decomposition[Feature Request]

[Recktenwald]

Is your feature request related to a problem? Please describe.
I noticed that the standard implementation was relatively slow. Here is a code example

open FSharp.Stats
open FSharp.Stats.Distributions
open FSharp.Stats.Fitting.LinearRegression
let normal = Continuous.normal 0.0 1.0
let n = 10000
let xs = [|for _ in [1..n] do 4.0 * normal.Sample()|]
let us = [|for _ in [1..n] do 36.0 * normal.Sample()|]
let yData = vector (Array.map2 (fun x u -> 3.0+2.0*x + u) xs us)
let xData = Matrix.ofJaggedArray (Array.zip xs us |> Array.map (fun (a,b) -> [|a;b|]))

#time
let olsCoeffs = OrdinaryLeastSquares.Linear.Multivariable.coefficients xData yData

olsCoeffs // Real: 00:00:10.390, CPU: 00:00:11.515, GC gen0: 143, gen1: 3, gen2: 2

#time
let X = Matrix.init (xData.NumRows) (xData.NumCols+1) (fun i j -> if j = 0 then 1.0 else xData.[i,j-1])
let XT = Matrix.transpose X
let lower = (XT * X) |> Algebra.LinearAlgebra.Cholesky
let gamma = Algebra.LinearAlgebra.SolveTriangularLinearSystem (Matrix.transpose lower) (XT * yData) false 
let beta = Algebra.LinearAlgebra.SolveTriangularLinearSystem  lower gamma true
beta // Real: 00:00:00.093, CPU: 00:00:00.000, GC gen0: 0, gen1: 0, gen2: 0

The result is already slightly unstable, the coefficients deviating further from the true value than in the standard implementation, but it is fast and useful for exploratory analysis.

Describe the solution you'd like
Add an implementation based on the Cholesky Decomposition.

Describe alternatives you've considered
Waiting longer for results?

Additional context
I have already done it and I'm only opening this issue because the contribution guidelines say to do it. Also as a side note the PR Template is leading to a 404 ( https://github.com/fslaborg/FSharp.Stats/blob/developer/PULL_REQUEST_TEMPLATE.md )