reduce allocations in ols and lts

Author: jbytecode

CHANGELOG.md

@@ -1,5 +1,6 @@
 # v0.11.8 (Upcoming Release)
 
+- Reduce allocations in OLS and LTS. Implements new olsf() and olsf!() functions.
 
 # v0.11.7 
 

src/LinRegOutliers.jl

@@ -63,7 +63,7 @@ import .DataSets: hills, softdrinkdelivery, animals
 # Ordinary least squares type and functions
 # for fast regression tasks in outlier detection algorithms
 include("ols.jl")
-import .OrdinaryLeastSquares: OLS, ols, wls, residuals, predict, coef
+import .OrdinaryLeastSquares: OLS, ols, wls, residuals, predict, coef, olsf, olsf!
 
 # Regression diagnostics
 include("diagnostics.jl")
@@ -250,7 +250,7 @@ export diagnose
 
 
 # Ordinary least squares
-export OLS, ols, wls, residuals, predict, coef
+export OLS, ols, wls, residuals, predict, coef, olsf, olsf!
 
 # Algorithms
 export hs93, hs93initialset, hs93basicsubset

src/lts.jl

@@ -4,7 +4,7 @@ export lts
 export iterateCSteps
 
 import ..Basis: RegressionSetting, @extractRegressionSetting, designMatrix, responseVector
-import ..OrdinaryLeastSquares: ols, residuals, coef
+import ..OrdinaryLeastSquares: residuals, coef, olsf!
 
 import Distributions: sample!, mean
 
@@ -50,9 +50,10 @@ function iterateCSteps(
     objective::Float64 = Inf64
     iter::Int = 0
     sortedresindices = Array{Int}(undef, n)
+    tempbetas = zeros(size(X, 2))
     while iter < maxiter
-        tempols = ols(view(X, subsetindices, :), view(y, subsetindices))
-        absres = abs.(y - X * coef(tempols))
+        olsf!(view(X, subsetindices, :), view(y, subsetindices), tempbetas)
+        absres = abs.(y - X * tempbetas)
         sortperm!(sortedresindices, absres)
         subsetindices = view(sortedresindices, 1:h)
         objective = sum(view(sort!(absres .^ 2.0), 1:h))

src/ols.jl

@@ -2,6 +2,7 @@ module OrdinaryLeastSquares
 
 
 export OLS, ols, wls, residuals, predict, coef
+export olsf, olsf!
 
 import LinearAlgebra: ColumnNorm, qr, Diagonal
 import ..Basis:
@@ -51,8 +52,14 @@ julia> reg.betas
 
 """
 ols(X::AbstractMatrix{Float64}, y::AbstractVector{Float64})::OLS = OLS(X, y, qr(X, ColumnNorm()) \ y)
-#ols(X::AbstractMatrix{Float64}, y::AbstractVector{Float64})::OLS = OLS(X, y, qr(X, Val(true)) \ y)
 
+function olsf(X::AbstractMatrix{Float64}, y::AbstractVector{Float64})::Vector{Float64}
+     return qr(X, ColumnNorm()) \ y
+end
+
+function olsf!(X::AbstractMatrix{Float64}, y::AbstractVector{Float64}, betas::AbstractVector{Float64})::Vector{Float64}
+     betas .= qr(X, ColumnNorm()) \ y
+end 
 
 function ols(setting::RegressionSetting)::OLS
     X, y = @extractRegressionSetting setting

test/testols.jl

@@ -23,6 +23,34 @@
         @test isapprox(result.betas[2], 5.04147826086957, atol = eps)
     end 
 
+    @testset "Ordinary Least Squares (olsf)" begin
+        tol = 0.0001
+        #  The model is exactly y = 5 + 5x
+        var1 = Float64[1, 2, 3, 4, 5]
+        X = hcat(ones(5), var1)
+        betas = [5.0, 5.0]
+        y = X * betas
+
+        # OLS
+        olsreg = olsf(X, y)
+        @test isapprox(olsreg, betas, atol = tol)
+    end
+
+    @testset "Ordinary Least Squares (olsf!)" begin
+        tol = 0.0001
+        #  The model is exactly y = 5 + 5x
+        var1 = Float64[1, 2, 3, 4, 5]
+        X = hcat(ones(5), var1)
+        betas = [5.0, 5.0]
+        y = X * betas
+
+        # OLS
+        tempbetas = zeros(size(X, 2))
+        olsf!(X, y, tempbetas)
+        @test isapprox(tempbetas, betas, atol = tol)
+    end
+
+
 
     @testset "Weighted Least Squares" begin
         tol = 0.0001