QR decomposition method in GLM.

docs/src/examples.md

@@ -86,6 +86,63 @@ julia> round.(vcov(ols); digits=5)
   0.38889  -0.16667
  -0.16667   0.08333
 ```
+By default, the `lm` method uses the Cholesky factorization which is known as fast but numerically unstable, especially for ill-conditioned design matrices. You can use the `method` keyword argument to apply a more stable QR factorization method.
+
+```jldoctest
+julia> data = DataFrame(X=[1,2,3], Y=[2,4,7]);
+
+julia> ols = lm(@formula(Y ~ X), data; method=:qr)
+LinearModel
+
+Y ~ 1 + X
+
+Coefficients:
+─────────────────────────────────────────────────────────────────────────
+                 Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
+─────────────────────────────────────────────────────────────────────────
+(Intercept)  -0.666667    0.62361   -1.07    0.4788   -8.59038    7.25704
+X             2.5         0.288675   8.66    0.0732   -1.16797    6.16797
+─────────────────────────────────────────────────────────────────────────
+```
+The following example shows that QR decomposition works better for an ill-conditioned design matrix. The linear model with the Cholesky decomposition method is unable to estimate parameters correctly whereas the linear model with the QR decomposition does.
+Note that, the condition number of the design matrix is quite high (≈ 3.52e7). Since the Cholesky method aggressively detect multicollinearity, if you ever encounter multicollinearity in any GLM model with Cholesky, it is worth trying the same model with QR decomposition. 
+
+```
+julia> X = [-0.4011512997627107 0.6368622664511552;
+            -0.0808472925693535 0.12835204623364604;
+            -0.16931095045225217 0.2687956795496601;
+            -0.4110745650568839 0.6526163576003452;
+            -0.4035951747670475 0.6407421349445884;
+            -0.4649907741370211 0.7382129928076485;
+            -0.15772708898883683 0.25040532268222715;
+            -0.38144358562952446 0.6055745630707645;
+            -0.1012787681395544 0.16078875117643368;
+            -0.2741403589052255 0.4352214984054432];
+
+julia> y = X * [5, 10];
+
+julia> md1 = lm(X, y; method=:qr, dropcollinear=false)
+LinearModel
+
+Coefficients:
+───────────────────────────────────────────────────────────────────
+    Coef.  Std. Error             t  Pr(>|t|)  Lower 95%  Upper 95%
+───────────────────────────────────────────────────────────────────
+x1    5.0  2.53054e-8  197586428.62    <1e-63        5.0        5.0
+x2   10.0  1.59395e-8  627370993.89    <1e-99       10.0       10.0
+───────────────────────────────────────────────────────────────────
+
+julia> md2 = lm(X, y; method=:cholesky, dropcollinear=false)
+LinearModel
+
+Coefficients:
+────────────────────────────────────────────────────────────────
+       Coef.  Std. Error       t  Pr(>|t|)  Lower 95%  Upper 95%
+────────────────────────────────────────────────────────────────
+x1   5.28865   0.103248    51.22    <1e-10    5.05056    5.52674
+x2  10.1818    0.0650348  156.56    <1e-14   10.0318    10.3318
+────────────────────────────────────────────────────────────────
+```
 
 ## Probit regression
 ```jldoctest

src/GLM.jl

@@ -91,13 +91,23 @@ module GLM
         pivoted_cholesky!(A; kwargs...) = cholesky!(A, RowMaximum(); kwargs...)
     end
 
+    @static if VERSION < v"1.7.0"
+        pivoted_qr!(A; kwargs...) = qr!(A, Val(true); kwargs...)
+    else
+        pivoted_qr!(A; kwargs...) = qr!(A, ColumnNorm(); kwargs...)
+    end
+
     const COMMON_FIT_KWARGS_DOCS = """
         - `dropcollinear::Bool=true`: Controls whether or not a model matrix
           less-than-full rank is accepted.
           If `true` (the default) the coefficient for redundant linearly dependent columns is
           `0.0` and all associated statistics are set to `NaN`.
           Typically from a set of linearly-dependent columns the last ones are identified as redundant
           (however, the exact selection of columns identified as redundant is not guaranteed).
+        - `method::Symbol=:cholesky`: Controls which decomposition method to use.
+          If `method=:cholesky` (the default), then the `Cholesky` decomposition method will be used.
+          If `method=:qr`, then the `QR` decomposition method (which is more stable
+          but slower) will be used.
         - `wts::Vector=similar(y,0)`: Prior frequency (a.k.a. case) weights of observations.
           Such weights are equivalent to repeating each observation a number of times equal
           to its weight. Do note that this interpretation gives equal point estimates but

src/glmfit.jl

@@ -286,7 +286,8 @@ function nulldeviance(m::GeneralizedLinearModel)
         X = fill(1.0, length(y), hasint ? 1 : 0)
         nullm = fit(GeneralizedLinearModel,
                     X, y, d, r.link; wts=wts, offset=offset,
-                    dropcollinear=isa(m.pp.chol, CholeskyPivoted),
+                    dropcollinear=ispivoted(m.pp),
+                    method=decomposition_method(m.pp),
                     maxiter=m.maxiter, minstepfac=m.minstepfac,
                     atol=m.atol, rtol=m.rtol)
         dev = deviance(nullm)
@@ -341,7 +342,8 @@ function nullloglikelihood(m::GeneralizedLinearModel)
         X = fill(1.0, length(y), hasint ? 1 : 0)
         nullm = fit(GeneralizedLinearModel,
                     X, y, d, r.link; wts=wts, offset=offset,
-                    dropcollinear=isa(m.pp.chol, CholeskyPivoted),
+                    dropcollinear=ispivoted(m.pp),
+                    method=decomposition_method(m.pp),
                     maxiter=m.maxiter, minstepfac=m.minstepfac,
                     atol=m.atol, rtol=m.rtol)
         ll = loglikelihood(nullm)
@@ -559,6 +561,7 @@ function fit(::Type{M},
     d::UnivariateDistribution,
     l::Link = canonicallink(d);
     dropcollinear::Bool = true,
+    method::Symbol = :cholesky,
     dofit::Union{Bool, Nothing} = nothing,
     wts::AbstractVector{<:Real}      = similar(y, 0),
     offset::AbstractVector{<:Real}   = similar(y, 0),
@@ -575,7 +578,15 @@ function fit(::Type{M},
     end
 
     rr = GlmResp(y, d, l, offset, wts)
-    res = M(rr, cholpred(X, dropcollinear), nothing, false)
+
+    if method === :cholesky
+        res = M(rr, cholpred(X, dropcollinear), nothing, false)
+    elseif method === :qr
+        res = M(rr, qrpred(X, dropcollinear), nothing, false)
+    else
+        throw(ArgumentError("The only supported values for keyword argument `method` are `:cholesky` and `:qr`."))
+    end
+
     return dofit ? fit!(res; fitargs...) : res
 end
 
@@ -594,6 +605,7 @@ function fit(::Type{M},
              offset::Union{AbstractVector, Nothing} = nothing,
              wts::Union{AbstractVector, Nothing} = nothing,
              dropcollinear::Bool = true,
+             method::Symbol = :cholesky,
              dofit::Union{Bool, Nothing} = nothing,
              contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}(),
              fitargs...) where {M<:AbstractGLM}
@@ -613,7 +625,15 @@ function fit(::Type{M},
     off = offset === nothing ? similar(y, 0) : offset
     wts = wts === nothing ? similar(y, 0) : wts
     rr = GlmResp(y, d, l, off, wts)
-    res = M(rr, cholpred(X, dropcollinear), f, false)
+
+    if method === :cholesky
+        res = M(rr, cholpred(X, dropcollinear), f, false)
+    elseif method === :qr
+        res = M(rr, qrpred(X, dropcollinear), f, false)
+    else
+        throw(ArgumentError("The only supported values for keyword argument `method` are `:cholesky` and `:qr`."))
+    end
+
     return dofit ? fit!(res; fitargs...) : res
 end
 

src/linpred.jl

@@ -35,45 +35,125 @@ end
 """
     DensePredQR
 
-A `LinPred` type with a dense, unpivoted QR decomposition of `X`
+A `LinPred` type with a dense QR decomposition of `X`
 
 # Members
 
 - `X`: Model matrix of size `n` × `p` with `n ≥ p`.  Should be full column rank.
 - `beta0`: base coefficient vector of length `p`
 - `delbeta`: increment to coefficient vector, also of length `p`
 - `scratchbeta`: scratch vector of length `p`, used in `linpred!` method
-- `qr`: a `QRCompactWY` object created from `X`, with optional row weights.
+- `qr`: either a `QRCompactWY` or `QRPivoted` object created from `X`, with optional row weights.
+- `scratchm1`: scratch Matrix{T} of the same size as `X`
 """
-mutable struct DensePredQR{T<:BlasReal} <: DensePred
+mutable struct DensePredQR{T<:BlasReal,Q<:Union{QRCompactWY, QRPivoted}} <: DensePred
     X::Matrix{T}                  # model matrix
     beta0::Vector{T}              # base coefficient vector
     delbeta::Vector{T}            # coefficient increment
     scratchbeta::Vector{T}
-    qr::QRCompactWY{T}
-    function DensePredQR{T}(X::Matrix{T}, beta0::Vector{T}) where T
+    qr::Q
+    scratchm1::Matrix{T}
+
+    function DensePredQR(X::AbstractMatrix, beta0::AbstractVector, pivot::Bool=false)
         n, p = size(X)
         length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
-        new{T}(X, beta0, zeros(T,p), zeros(T,p), qr(X))
+        T = typeof(float(zero(eltype(X))))
+        Q = pivot ? QRPivoted : QRCompactWY
+        fX = float(X)
+        cfX = fX === X ? copy(fX) : fX
+        F = pivot ? pivoted_qr!(copy(cfX)) : qr!(cfX)
+        new{T,Q}(Matrix{T}(X),
+            Vector{T}(beta0),
+            zeros(T, p),
+            zeros(T, p),
+            F,
+            similar(X, T))
     end
-    function DensePredQR{T}(X::Matrix{T}) where T
+    function DensePredQR(X::AbstractMatrix, pivot::Bool=false)
         n, p = size(X)
-        new{T}(X, zeros(T, p), zeros(T,p), zeros(T,p), qr(X))
+        T = typeof(float(zero(eltype(X))))
+        Q = pivot ? QRPivoted : QRCompactWY
+        F = pivot ? pivoted_qr!(float(copy(X))) : qr(float(X))
+        new{T,Q}(Matrix{T}(X),
+            zeros(T, p),
+            zeros(T, p),
+            zeros(T, p),
+            F,
+            similar(X, T))
     end
 end
-DensePredQR(X::Matrix, beta0::Vector) = DensePredQR{eltype(X)}(X, beta0)
-DensePredQR(X::Matrix{T}) where T = DensePredQR{T}(X, zeros(T, size(X,2)))
-convert(::Type{DensePredQR{T}}, X::Matrix{T}) where {T} = DensePredQR{T}(X, zeros(T, size(X, 2)))
-
 """
     delbeta!(p::LinPred, r::Vector)
 
 Evaluate and return `p.delbeta` the increment to the coefficient vector from residual `r`
 """
 function delbeta! end
 
-function delbeta!(p::DensePredQR{T}, r::Vector{T}) where T<:BlasReal
+function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}) where T<:BlasReal
+    rnk = rank(p.qr.R)
+    rnk == length(p.delbeta) || throw(RankDeficientException(rnk))
     p.delbeta = p.qr\r
+    mul!(p.scratchm1, Diagonal(ones(size(r))), p.X)
+    return p
+end
+
+function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal
+    rnk = rank(p.qr.R)
+    rnk == length(p.delbeta) || throw(RankDeficientException(rnk))
+    X = p.X
+    W = Diagonal(wt)
+    sqrtW = Diagonal(sqrt.(wt)) 
+    mul!(p.scratchm1, sqrtW, X)
+    mul!(p.delbeta, X'W, r)
+    qnr = qr(p.scratchm1)
+    Rinv = inv(qnr.R)
+    p.delbeta = Rinv * Rinv' * p.delbeta
+    return p
+end
+
+function delbeta!(p::DensePredQR{T,<:QRPivoted}, r::Vector{T}) where T<:BlasReal
+    rnk = rank(p.qr.R)
+    if rnk == length(p.delbeta)
+        p.delbeta = p.qr\r
+    else
+        R = @view p.qr.R[:, 1:rnk] 
+        Q = @view p.qr.Q[:, 1:size(R, 1)]
+        piv = p.qr.p
+        p.delbeta = zeros(size(p.delbeta))
+        p.delbeta[1:rnk] = R \ Q'r
+        invpermute!(p.delbeta, piv)
+    end
+    mul!(p.scratchm1, Diagonal(ones(size(r))), p.X)
+    return p
+end
+
+function delbeta!(p::DensePredQR{T,<:QRPivoted}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal
+    rnk = rank(p.qr.R)
+    X = p.X
+    W = Diagonal(wt)
+    sqrtW = Diagonal(sqrt.(wt))
+    delbeta = p.delbeta
+    scratchm2 = similar(X, T)
+    mul!(p.scratchm1, sqrtW, X)
+    mul!(scratchm2, W, X)
+    mul!(delbeta, transpose(scratchm2), r)
+
+    if rnk == length(p.delbeta)
+        qnr = qr(p.scratchm1)
+        Rinv = inv(qnr.R)
+        p.delbeta = Rinv * Rinv' * delbeta
+    else
+        qnr = pivoted_qr!(copy(p.scratchm1))
+        R = @view qnr.R[1:rnk, 1:rnk]
+        Rinv = inv(R)
+        piv = qnr.p
+        permute!(delbeta, piv)
+        for k=(rnk+1):length(delbeta)
+            delbeta[k] = -zero(T)
+        end
+        p.delbeta[1:rnk] = Rinv * Rinv' * view(delbeta, 1:rnk)
+        invpermute!(delbeta, piv)
+    end
     return p
 end
 
@@ -115,6 +195,7 @@ function DensePredChol(X::AbstractMatrix, pivot::Bool)
 end
 
 cholpred(X::AbstractMatrix, pivot::Bool=false) = DensePredChol(X, pivot)
+qrpred(X::AbstractMatrix, pivot::Bool=false) = DensePredQR(X, pivot)
 
 cholfactors(c::Union{Cholesky,CholeskyPivoted}) = c.factors
 cholesky!(p::DensePredChol{T}) where {T<:FP} = p.chol
@@ -124,7 +205,6 @@ function cholesky(p::DensePredChol{T}) where T<:FP
     c = p.chol
     Cholesky(copy(cholfactors(c)), c.uplo, c.info)
 end
-cholesky!(p::DensePredQR{T}) where {T<:FP} = Cholesky{T,typeof(p.X)}(p.qr.R, 'U', 0)
 
 function delbeta!(p::DensePredChol{T,<:Cholesky}, r::Vector{T}) where T<:BlasReal
     ldiv!(p.chol, mul!(p.delbeta, transpose(p.X), r))
@@ -227,7 +307,33 @@ end
 LinearAlgebra.cholesky(p::SparsePredChol{T}) where {T} = copy(p.chol)
 LinearAlgebra.cholesky!(p::SparsePredChol{T}) where {T} = p.chol
 
+function invqr(x::DensePredQR{T,<: QRCompactWY}) where T
+    Q,R = qr(x.scratchm1)
+    Rinv = inv(R)
+    Rinv*Rinv'
+end
+
+function invqr(x::DensePredQR{T,<: QRPivoted}) where T
+    Q,R,pv = pivoted_qr!(copy(x.scratchm1))
+    rnk = rank(R)
+    p = length(x.delbeta)
+    if rnk == p
+        Rinv = inv(R)
+        xinv = Rinv*Rinv'
+        ipiv = invperm(pv)
+        return xinv[ipiv, ipiv]
+    else
+        Rsub = R[1:rnk, 1:rnk]
+        RsubInv = inv(Rsub)
+        xinv = fill(convert(T, NaN), (p,p))
+        xinv[1:rnk, 1:rnk] = RsubInv*RsubInv'
+        ipiv = invperm(pv)
+        return xinv[ipiv, ipiv]
+    end
+end
+
 invchol(x::DensePred) = inv(cholesky!(x))
+
 function invchol(x::DensePredChol{T,<: CholeskyPivoted}) where T
     ch = x.chol
     rnk = rank(ch)
@@ -242,8 +348,14 @@ function invchol(x::DensePredChol{T,<: CholeskyPivoted}) where T
     ipiv = invperm(ch.p)
     res[ipiv, ipiv]
 end
+
 invchol(x::SparsePredChol) = cholesky!(x) \ Matrix{Float64}(I, size(x.X, 2), size(x.X, 2))
-vcov(x::LinPredModel) = rmul!(invchol(x.pp), dispersion(x, true))
+
+inverse(x::DensePred) = invchol(x)
+inverse(x::DensePredQR) = invqr(x)
+inverse(x::SparsePredChol) = invchol(x)
+
+vcov(x::LinPredModel) = rmul!(inverse(x.pp), dispersion(x, true))
 
 function cor(x::LinPredModel)
     Σ = vcov(x)
@@ -308,7 +420,15 @@ dof_residual(obj::LinPredModel) = nobs(obj) - dof(obj) + 1
 hasintercept(m::LinPredModel) = any(i -> all(==(1), view(m.pp.X , :, i)), 1:size(m.pp.X, 2))
 
 linpred_rank(x::LinPred) = length(x.beta0)
-linpred_rank(x::DensePredChol{<:Any, <:CholeskyPivoted}) = x.chol.rank
+linpred_rank(x::DensePredChol{<:Any, <:CholeskyPivoted}) = rank(x.chol)
+linpred_rank(x::DensePredChol{<:Any, <:Cholesky}) = rank(x.chol.U)
+linpred_rank(x::DensePredQR{<:Any,<:QRPivoted}) = rank(x.qr.R)
+
+ispivoted(x::LinPred) = false
+ispivoted(x::DensePredChol{<:Any, <:CholeskyPivoted}) = true
+ispivoted(x::DensePredQR{<:Any,<:QRPivoted}) = true
+
+decomposition_method(x::LinPred) = isa(x, DensePredQR) ? :qr : :cholesky
 
 _coltype(::ContinuousTerm{T}) where {T} = T