I was inspecting GLM.jl IRLS implementation, but haven’t figured out from the code how many matrix decompositions it performs.

For any given, `A::AbstractMatrix{<:Real}`

, `b::AbstractVector{<:Real}`

, and `wts::AbstractWeights`

,

for solving the linear system: Ax = b with weights w, x = (A^{\top}WA)^{-1}A^{\top}Wb.

Assuming I compute the `F::QRCompactWY = qrfact(sqrt.(w) .* A)`

and solve `x = F \ (sqrt.(w) .* y)`

, for new weights w_{2}, how can I solve the new system relying on my already computed `F::QRCompactWY`

? Do I need to compute a new factorization every iteration? This study suggests the best would be a Q-less decomposition.