Benchmarking etc (#29)

* decreasing allocs via temp

* using morel legible refs

* reviews 1

* adding tests for scratchspace

* add scratch for robust

* tests ok

* adding a bunch of tests and completing the PR

* additional fixes

* update bench

* fixing all conflicts...

Author: tlienart

src/glr/d_l2loss.jl

@@ -25,7 +25,7 @@ function Hv!(glr::GLR{L2Loss,<:L2R}, X, y)
             Hvₐ = view(Hv, a)
             vₐ  = view(v,  a)
             Xt1 = view(SCRATCH_P[], a)
-            copyto!(Xt1, sum(X, dims=1))  # -- X'1
+            copyto!(Xt1, sum(X, dims=1))  # -- X'1 (note: sum will allocate)
             vₑ  = v[end]
             # update for the first p rows   -- (X'X + λI)v[1:p] + (X'1)v[end]
             Xvₐ = SCRATCH_N[]
@@ -60,10 +60,8 @@ function smooth_fg!(glr::GLR{L2Loss,<:ENR}, X, y)
     λ = getscale_l2(glr.penalty)
     (g, θ) -> begin
         # cache contains the residuals (Xθ-y)
-        Xθ = SCRATCH_N[]
-        apply_X!(Xθ, X, θ)
-        r   = SCRATCH_N[]
-        r .-= y             # -- r = Xθ-y
+        r = SCRATCH_N[]
+        get_residuals!(r, X, θ, y) # -- r = Xθ-y
         apply_Xt!(g, X, r)
         g .+= λ .* θ
         return glr.loss(r) + get_l2(glr.penalty)(θ)

src/glr/d_robust.jl

@@ -38,17 +38,20 @@ function fgh!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y) where ρ <: RobustRho1P{δ}
     ϕ_ = ϕ(ρ)
     if glr.fit_intercept
         (f, g, H, θ) -> begin
-            r = _get_residuals(X, θ, y)
-            w = _get_w(r, δ)
+            r  = SCRATCH_N[]
+            get_residuals!(r, X, θ, y)
+            w  = SCRATCH_N2[]
+            w .= convert.(Float64, abs.(r) .<= δ)
             # gradient via ψ function
             g === nothing || begin
-                ψr = _get_ψr(r, w, ψ_)
+                ψr  = SCRATCH_N3[]
+                ψr .= ψ_.(r, w)
                 apply_Xt!(g, X, ψr)
                 g .+= λ .* θ
             end
             # Hessian via ϕ functiono
             H === nothing || begin
-                # Hessian allocates a ton anyway so use of scratch is a bit pointless
+                # NOTE: Hessian allocates a ton anyway so use of scratch is a bit pointless
                 ϕr = ϕ_.(r, w)
                 ΛX = ϕr .* X
                 mul!(view(H, 1:p, 1:p), X', ΛX)
@@ -64,11 +67,14 @@ function fgh!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y) where ρ <: RobustRho1P{δ}
         end
     else
         (f, g, H, θ) -> begin
-            r = _get_residuals(X, θ, y)
-            w = _get_w(r, δ)
+            r = SCRATCH_N[]
+            get_residuals!(r, X, θ, y)
+            w = SCRATCH_N2[]
+            w .= convert.(Float64, abs.(r) .<= δ)
             # gradient via ψ function
             g === nothing || begin
-                ψr = _get_ψr(r, w, ψ_)
+                ψr  = SCRATCH_N3[]
+                ψr .= ψ_.(r, w)
                 apply_Xt!(g, X, ψr)
                 g .+= λ .* θ
             end
@@ -87,8 +93,10 @@ function Hv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y) where ρ <: RobustRho1P{δ} w
     # see d_logistic.jl for more comments on this (similar procedure)
     if glr.fit_intercept
         (Hv, θ, v) -> begin
-            r  = _get_residuals(X, θ, y)
-            w  = _get_w(r, δ)
+            r  = SCRATCH_N[]
+            get_residuals!(r, X, θ, y)
+            w  = SCRATCH_N2[]
+            w .= convert.(Float64, abs.(r) .<= δ)
             w .= ϕ_.(r, w)
             # views on first p rows (intercept row treated after)
             a    = 1:p
@@ -108,8 +116,10 @@ function Hv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y) where ρ <: RobustRho1P{δ} w
         end
     else
         (Hv, θ, v) -> begin
-            r  = _get_residuals(X, θ, y)
-            w  = _get_w(r, δ)
+            r  = SCRATCH_N[]
+            get_residuals!(r, X, θ, y)
+            w  = SCRATCH_N2[]
+            w .= convert.(Float64, abs.(r) .<= δ)
             w .= ϕ_.(r, w)
             t  = SCRATCH_N3[]
             apply_X!(t, X, v)
@@ -120,6 +130,7 @@ function Hv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y) where ρ <: RobustRho1P{δ} w
     end
 end
 
+
 # For IWLS
 function Mv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y;
              threshold=1e-6) where ρ <: RobustRho1P{δ} where δ
@@ -130,24 +141,38 @@ function Mv!(glr::GLR{RobustLoss{ρ},<:L2R}, X, y;
     # which we solve via an iterative method so, one θ
     # gives one way of applying the relevant matrix (X'ΛX+λI)
     (ωr, θ) -> begin
-        r = _get_residuals(X, θ, y)
-        w = _get_w(r, δ)
+        r   = SCRATCH_N[]
+        get_residuals!(r, X, θ, y)
+        w   = SCRATCH_N2[]
+        w  .= convert.(Float64, abs.(r) .<= δ)
         # ω = ψ(r)/r ; weighing factor for IWLS
         ωr .= ω_.(r, w)
         # function defining the application of (X'ΛX + λI)
         if glr.fit_intercept
             (Mv, v) -> begin
-                a    = 1:p
-                vₐ   = view(v, a)
-                Mvₐ  = view(Mv, a)
-                XtW1 = vec(sum(ωr .* X, dims=1))
-                vₑ   = v[end]
-                mul!(Mvₐ, X', ωr .* (X * vₐ))
+                a     = 1:p
+                vₐ    = view(v, a)
+                Mvₐ   = view(Mv, a)
+                XtW1  = view(SCRATCH_P[], a)
+                @inbounds for j in a
+                    XtW1[j] = dot(ωr, view(X, :, j))
+                end
+                vₑ = v[end]
+                t  = SCRATCH_N[]
+                apply_X!(t, X, vₐ)
+                t .*= ωr
+                mul!(Mvₐ, X', t)
                 Mvₐ .+= λ .* vₐ .+ XtW1 .* vₑ
                 Mv[end] = dot(XtW1, vₐ) + (sum(ωr)+λ) * vₑ
             end
         else
-            (Mv, v) -> (mul!(Mv, X', ωr .* (X * v));  Mv .+= λ .* v)
+            (Mv, v) -> begin
+                t  = SCRATCH_N[]
+                apply_X!(t, X, v)
+                t .*= ωr
+                mul!(Mv, X', t)
+                Mv .+= λ .* v
+            end
         end
     end
 end
@@ -159,9 +184,12 @@ function smooth_fg!(glr::GLR{RobustLoss{ρ},<:ENR}, X, y) where ρ <: RobustRho1
     p  = size(X, 2)
     ψ_ = ψ(ρ)
     (g, θ) -> begin
-        r  = _get_residuals(X, θ, y)
-        w  = _get_w(r, δ)
-        ψr = _get_ψr(r, w, ψ_)
+        r   = SCRATCH_N[]
+        get_residuals!(r, X, θ, y)
+        w   = SCRATCH_N2[]
+        w  .= convert.(Float64, abs.(r) .<= δ)
+        ψr  = SCRATCH_N3[]
+        ψr .= ψ_.(r, w)
         apply_Xt!(g, X, ψr)
         g .+= λ .* θ
         return glr.loss(r) + get_l2(glr.penalty)(θ)

src/glr/utils.jl

@@ -40,3 +40,14 @@ $SIGNATURES
 Return a model corresponding to the smooth part of the objective.
 """
 get_smooth(glr::GLR) = (o = smooth_objective(glr); GLR(o.loss, o.penalty, glr.fit_intercept))
+
+
+"""
+$SIGNATURES
+
+Helper function to compute the residuals.
+"""
+function get_residuals!(r, X, θ, y)
+    apply_X!(r, X, θ)
+    r .-= y
+end

src/utils.jl

@@ -57,7 +57,7 @@ function apply_X!(Xθ, X, θ, c=1)
 		if length(θ) == p
 			mul!(Xθ, X, θ)
 		else
-			mul!(Xθ, X, θ[1:p])
+			mul!(Xθ, X, view(θ, 1:p))
 			Xθ .+= θ[end]
 		end
 	else

test/benchmarks/elementary_functions.jl

@@ -0,0 +1,50 @@
+using MLJLinearModels
+using BenchmarkTools, Random, LinearAlgebra
+DO_COMPARISONS = false; include("../testutils.jl")
+
+n, p = 50_000, 500
+((X, y, θ), (X_, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)
+
+# ======================== #
+# ELEMENTARIES
+#
+# > apply_X!        ✅ Sept 13, 2019
+# > apply_Xt!       ✅ Sept 13, 2019
+# > get_residuals!  ✅ Sept 13, 2019
+# ======================== #
+
+# ------------
+# No intercept
+
+R.allocate(n, p)
+
+Xθ = similar(y)
+Xtθ = similar(θ)
+r = R.SCRATCH_N[]
+
+# Sept 13, 2019 :: 7.83 ms (0 allocations: 0 bytes)
+@btime R.apply_X!($Xθ, $X, $θ);
+
+# Sept 13, 2019 :: 7.94 ms (0 allocations: 0 bytes) [pretty much only apply_X!]
+@btime R.get_residuals!($r, $X, $θ, $y);
+
+# Sept 13, 2019 :: 7.9 ms (0 allocations: 0 bytes)
+@btime R.apply_Xt!($Xtθ, $X, $y);
+
+# --------------
+# With intercept
+
+R.allocate(n, p+1)
+
+Xθ = similar(y)
+Xtθ = similar(θ1)
+r = R.SCRATCH_N[]
+
+# Sept 13, 2019 :: 7.949 ms (1 allocation: 48 bytes) -- alloc for the view
+@btime R.apply_X!($Xθ, $X, $θ1);
+
+# Sept 13, 2019 :: 7.955 ms (1 allocation: 48 bytes) -- mostly apply_X!
+@btime R.get_residuals!($r, $X, $θ1, $y);
+
+# Sept 13, 2019 :: 7.883 ms (1 allocation: 48 bytes)
+@btime R.apply_Xt!($Xtθ, $X, $y);

test/benchmarks/ridge-lasso.jl

@@ -1,15 +1,74 @@
-# WIP WIP !
-
 using MLJLinearModels
 using BenchmarkTools, Random, LinearAlgebra
 DO_COMPARISONS = false; include("../testutils.jl")
 
 n, p = 50_000, 500
+((X, y, θ), (X_, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)
+
+# =============== #
+# RIDGE FUNCTIONS #
+#
+# Hv!           ✅
+# =============== #
+
+# No fit_intercept
+R.allocate(n, p)
+ridge = RidgeRegression(0.5; fit_intercept=false)
+Hv! = R.Hv!(ridge, X, y)
+v   = randn(p)
+Hv  = similar(v)
+
+# Sept 13, 2019 :: 15.300 ms (0 allocations: 0 bytes)
+@btime Hv!($Hv, $v)
+
+# With fit_intercept
+R.allocate(n, p+1)
+ridge = RidgeRegression(λ)
+Hv! = R.Hv!(ridge, X, y)
+v  = randn(p+1)
+Hv = similar(v)
+
+# Sept 13, 2019 :: 26.246 ms (5 allocations: 4.22 KiB)
+# (decrease in perf due to views, and need a sum on an array)
+@btime Hv!($Hv, $v)
+
+# =============== #
+# LASSO FUNCTIONS #
+#
+# smooth_fg!
+# =============== #
+
+# No fit_intercept
+R.allocate(n, p)
+
+lasso = LassoRegression(0.5, fit_intercept=false)
+smooth_fg! = R.smooth_fg!(lasso, X, y)
+v = randn(p)
+g = similar(v)
+
+# Sept 13, 2019 :: 15.508 ms (3 allocations: 390.72 KiB)
+@btime smooth_fg!($g, $v);
+
+@btime R.get_residuals!($X, $θ, $y)
+
+r = R.get_residuals!(X, θ, y)
+
+@btime R.apply_Xt!(g, X, r)
+
+@btime (g .+= 0.5 .* θ)
+
+
+
+
+
+# =======
+=======
 ((X, _, _), (X_, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.1)
 
 λ  = 50
 lr = LassoRegression(λ)
 J  = objective(lr, X, y1)
+>>>>>>> master
 
 fista = FISTA()
 ista = ISTA()