Skip to content

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

@Recktenwald

Description

@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 )

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions