Allow for sparse arrays and views (#446)

JuliaStats · Jan 24, 2022 · fa74d14 · fa74d14 · andreasnoack · Jan 24, 2022
1 parent be96833
commit fa74d14
Show file tree

Hide file tree

Showing 5 changed files with 91 additions and 15 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "GLM"
 uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
-version = "1.6.0"
+version = "1.6.1"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

diff --git a/src/glmfit.jl b/src/glmfit.jl
@@ -49,11 +49,15 @@ function GlmResp(y::V, d::D, l::L, η::V, μ::V, off::V, wts::V) where {V<:FPVec
     return GlmResp{V,D,L}(y, d, similar(y), η, μ, off, wts, similar(y), similar(y))
 end
 
-function GlmResp(y::V, d::D, l::L, off::V, wts::V) where {V<:FPVector,D,L}
-    η   = similar(y)
-    μ   = similar(y)
-    r   = GlmResp(y, d, l, η, μ, off, wts)
-    initialeta!(r.eta, d, l, y, wts, off)
+function GlmResp(y::FPVector, d::Distribution, l::Link, off::FPVector, wts::FPVector)
+    # Instead of convert(Vector{Float64}, y) to be more ForwardDiff friendly
+    _y   = convert(Vector{float(eltype(y))}, y)
+    _off = convert(Vector{float(eltype(off))}, off)
+    _wts = convert(Vector{float(eltype(wts))}, wts)
+    η    = similar(_y)
+    μ    = similar(_y)
+    r    = GlmResp(_y, d, l, η, μ, _off, _wts)
+    initialeta!(r.eta, d, l, _y, _wts, _off)
     updateμ!(r, r.eta)
     return r
 end
@@ -465,14 +469,14 @@ Fit a generalized linear model to data.
 $FIT_GLM_DOC
 """
 function fit(::Type{M},
-    X::Union{Matrix{T},SparseMatrixCSC{T}},
+    X::AbstractMatrix{<:FP},
     y::AbstractVector{<:Real},
     d::UnivariateDistribution,
     l::Link = canonicallink(d);
     dofit::Bool = true,
     wts::AbstractVector{<:Real}      = similar(y, 0),
     offset::AbstractVector{<:Real}   = similar(y, 0),
-    fitargs...) where {M<:AbstractGLM,T<:FP}
+    fitargs...) where {M<:AbstractGLM}
 
     # Check that X and y have the same number of observations
     if size(X, 1) != size(y, 1)
@@ -485,7 +489,7 @@ function fit(::Type{M},
 end
 
 fit(::Type{M},
-    X::Union{Matrix,SparseMatrixCSC},
+    X::AbstractMatrix,
     y::AbstractVector,
     d::UnivariateDistribution,
     l::Link=canonicallink(d); kwargs...) where {M<:AbstractGLM} =

diff --git a/src/linpred.jl b/src/linpred.jl
@@ -185,7 +185,7 @@ function SparsePredChol(X::SparseMatrixCSC{T}) where T
         similar(X))
 end
 
-cholpred(X::SparseMatrixCSC) = SparsePredChol(X)
+cholpred(X::SparseMatrixCSC, pivot::Bool=false) = SparsePredChol(X)
 
 function delbeta!(p::SparsePredChol{T}, r::Vector{T}, wt::Vector{T}) where T
     scr = mul!(p.scratch, Diagonal(wt), p.X)
@@ -194,6 +194,13 @@ function delbeta!(p::SparsePredChol{T}, r::Vector{T}, wt::Vector{T}) where T
     p.delbeta = c \ mul!(p.delbeta, adjoint(scr), r)
 end
 
+function delbeta!(p::SparsePredChol{T}, r::Vector{T}) where T
+    scr = p.scratch = p.X
+    XtWX = p.Xt*scr
+    c = p.chol = cholesky(Symmetric{eltype(XtWX),typeof(XtWX)}(XtWX, 'L'))
+    p.delbeta = c \ mul!(p.delbeta, adjoint(scr), r)
+end
+
 LinearAlgebra.cholesky(p::SparsePredChol{T}) where {T} = copy(p.chol)
 LinearAlgebra.cholesky!(p::SparsePredChol{T}) where {T} = p.chol
 

diff --git a/src/lm.jl b/src/lm.jl
@@ -28,11 +28,16 @@ mutable struct LmResp{V<:FPVector} <: ModResp  # response in a linear model
     end
 end
 
-LmResp(y::FPVector, wts::FPVector=similar(y, 0)) = 
-    LmResp{typeof(y)}(fill!(similar(y), 0), similar(y, 0), wts, y)
-
-LmResp(y::AbstractVector{<:Real}, wts::AbstractVector{<:Real}=similar(y, 0)) = 
-    LmResp(float(y), float(wts))
+function LmResp(y::AbstractVector{<:Real}, wts::Union{Nothing,AbstractVector{<:Real}}=nothing)
+    # Instead of convert(Vector{Float64}, y) to be more ForwardDiff friendly
+    _y = convert(Vector{float(eltype(y))}, y)
+    _wts = if wts === nothing
+        similar(_y, 0)
+    else
+        convert(Vector{float(eltype(wts))}, wts)
+    end
+    return LmResp{typeof(_y)}(zero(_y), zero(_y), _wts, _y)
+end
 
 function updateμ!(r::LmResp{V}, linPr::V) where V<:FPVector
     n = length(linPr)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -502,6 +502,23 @@ end
     @test isapprox(vcov(gmsparse), vcov(gmdense))
 end
 
+@testset "Sparse LM" begin
+    rng = StableRNG(1)
+    X = sprand(rng, 1000, 10, 0.01)
+    β = randn(rng, 10)
+    y = Bool[rand(rng) < logistic(x) for x in X * β]
+    gmsparsev = [fit(LinearModel, X, y),
+                 fit(LinearModel, X, sparse(y)),
+                 fit(LinearModel, Matrix(X), sparse(y))]
+    gmdense = fit(LinearModel, Matrix(X), y)
+
+    for gmsparse in gmsparsev
+        @test isapprox(deviance(gmsparse), deviance(gmdense))
+        @test isapprox(coef(gmsparse), coef(gmdense))
+        @test isapprox(vcov(gmsparse), vcov(gmdense))
+    end
+end
+
 @testset "Predict" begin
     rng = StableRNG(123)
     X = rand(rng, 10, 2)
@@ -969,3 +986,46 @@ end
     secondcolinterceptmod = glm([randn(rng, 5) ones(5)], ones(5), Binomial(), LogitLink())
     @test hasintercept(secondcolinterceptmod)
 end
+
+@testset "Issue #444. Views" begin
+    X = randn(10, 2)
+    y = X*ones(2) + randn(10)
+    @test coef(glm(X, y, Normal(), IdentityLink())) ==
+        coef(glm(view(X, 1:10, :), view(y, 1:10), Normal(), IdentityLink()))
+
+    x, y, w = rand(100, 2), rand(100), rand(100)
+    lm1 = lm(x, y)
+    lm2 = lm(x, view(y, :))
+    lm3 = lm(view(x, :, :), y)
+    lm4 = lm(view(x, :, :), view(y, :))
+    @test coef(lm1) == coef(lm2) == coef(lm3) == coef(lm4)
+
+    lm5 = lm(x, y, wts=w)
+    lm6 = lm(x, view(y, :), wts=w)
+    lm7 = lm(view(x, :, :), y, wts=w)
+    lm8 = lm(view(x, :, :), view(y, :), wts=w)
+    lm9 = lm(x, y, wts=view(w, :))
+    lm10 = lm(x, view(y, :), wts=view(w, :))
+    lm11 = lm(view(x, :, :), y, wts=view(w, :))
+    lm12 = lm(view(x, :, :), view(y, :), wts=view(w, :))
+    @test coef(lm5) == coef(lm6) == coef(lm7) == coef(lm8) == coef(lm9) == coef(lm10) ==
+        coef(lm11) == coef(lm12)
+
+    x, y, w = rand(100, 2), rand(Bool, 100), rand(100)
+    glm1 = glm(x, y, Binomial())
+    glm2 = glm(x, view(y, :), Binomial())
+    glm3 = glm(view(x, :, :), y, Binomial())
+    glm4 = glm(view(x, :, :), view(y, :), Binomial())
+    @test coef(glm1) == coef(glm2) == coef(glm3) == coef(glm4)
+
+    glm5 = glm(x, y, Binomial(), wts=w)
+    glm6 = glm(x, view(y, :), Binomial(), wts=w)
+    glm7 = glm(view(x, :, :), y, Binomial(), wts=w)
+    glm8 = glm(view(x, :, :), view(y, :), Binomial(), wts=w)
+    glm9 = glm(x, y, Binomial(), wts=view(w, :))
+    glm10 = glm(x, view(y, :), Binomial(), wts=view(w, :))
+    glm11 = glm(view(x, :, :), y, Binomial(), wts=view(w, :))
+    glm12 = glm(view(x, :, :), view(y, :), Binomial(), wts=view(w, :))
+    @test coef(glm5) == coef(glm6) == coef(glm7) == coef(glm8) == coef(glm9) == coef(glm10) ==
+        coef(glm11) == coef(glm12)
+end