Skip to content

Rank calculation can vary during fitting of a GeneralizedLinearModel #601

Description

@andreasnoack

The idea in

GLM.jl/src/linpred.jl

Lines 204 to 222 in fb6c7a9

rnk = rank(p.chol)
if rnk == length(delbeta)
cf = cholfactors(p.chol)
cf .= p.scratchm2[piv, piv]
cholesky!(Hermitian(cf, Symbol(p.chol.uplo)))
ldiv!(p.chol, delbeta)
else
permute!(delbeta, piv)
for k=(rnk+1):length(delbeta)
delbeta[k] = -zero(T)
end
# shift full rank column to 1:rank
cf = cholfactors(p.chol)
cf .= p.scratchm2[piv, piv]
cholesky!(Hermitian(view(cf, 1:rnk, 1:rnk), Symbol(p.chol.uplo)))
ldiv!(Cholesky(view(cf, 1:rnk, 1:rnk), Symbol(p.chol.uplo), p.chol.info),
view(delbeta, 1:rnk))
invpermute!(delbeta, piv)
end
is to determine that rank and permutation vector for the pivoted Cholesky factorization before starting the IRLS algorithm and then apply the same permutations and rank during all subsequent iterations. However, as seen in the first line, the rank is recomputed at each iteration. This is inconsistent with fixing the rank related fields and could lead to invalid results. It might not be easy to produce are failing example, though.

For completeness, I should should mention that it isn't feasible to recompute the pivots and rank information during the optimization as it would violate the IRLS algorithm to zero out different coordinates at different iterations.

This came up when working on #600

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