-
Notifications
You must be signed in to change notification settings - Fork 114
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fix QR based fitting #559
base: master
Are you sure you want to change the base?
Fix QR based fitting #559
Conversation
Avoid erroring out for low rank design matrices when dropcollinear=false. Avoid unnecessary triangular solves. Avoid indexing in the Q. Avoid slicing R matrix in a way that triggers a minimum norm solution. Remove unnecessary temporaries.
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #559 +/- ##
==========================================
- Coverage 90.04% 89.87% -0.18%
==========================================
Files 8 8
Lines 1125 1106 -19
==========================================
- Hits 1013 994 -19
Misses 112 112 ☔ View full report in Codecov by Sentry. |
fill!(p.delbeta, 0) | ||
p.delbeta[1:rnk] = R \ view(p.qr.Q'r, 1:rnk) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe one could save a few allocations here?
fill!(p.delbeta, 0) | |
p.delbeta[1:rnk] = R \ view(p.qr.Q'r, 1:rnk) | |
mul!(view(p.delbeta, 1:rnk), view(p.qr.Q, :, 1:rnk)', r) | |
ldiv!(R, view(p.delbeta, 1:rnk)) | |
fill!(@view(p.delbeta[(rnk + 1):end]), 0) |
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is also unused?
W = Diagonal(wt) |
@@ -63,23 +63,19 @@ function delbeta! end | |||
|
|||
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are scratch-spaces only used inside of a function but not carried around? Or is there another reason for why this is not needed anymore?
qnr = qr(p.scratchm1) | ||
Rinv = inv(qnr.R) | ||
p.delbeta = Rinv * Rinv' * p.delbeta | ||
ỹ = sqrtW * r |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think one could avoid this allocation by keeping sqrt.(wt)
around separately,
sqrt_wt = sqrt.(wt)
mul!(p.scratchm1, Diagonal(sqrt_wt), X)
ỹ = (sqrt_wt .*= r)
@@ -88,44 +84,32 @@ function delbeta!(p::DensePredQR{T,<:QRPivoted}, r::Vector{T}) where T<:BlasReal | |||
if rnk == length(p.delbeta) | |||
p.delbeta = p.qr\r |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe
p.delbeta = p.qr\r | |
ldiv!(p.delbeta, p.qr, r) |
?
end | ||
p.delbeta[1:rnk] = Rinv * Rinv' * view(delbeta, 1:rnk) | ||
invpermute!(delbeta, piv) | ||
r̃ = sqrtW * r |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems the same approach as above could applied here?
invpermute!(delbeta, piv) | ||
r̃ = sqrtW * r | ||
|
||
p.qr = pivoted_qr!(copy(p.scratchm1)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is scratchm1
reused in another function?
R = UpperTriangular(view(parent(p.qr.R), 1:rnk, 1:rnk)) | ||
permute!(p.delbeta, p.qr.p) | ||
for k = (rnk + 1):length(p.delbeta) | ||
p.delbeta[k] = -zero(T) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems that's copied from the current implementation - but why -zero(T)
and not zero(T)
?
end | ||
p.delbeta[1:rnk] = R \ (p.qr.Q'*r̃)[1:rnk] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe
p.delbeta[1:rnk] = R \ (p.qr.Q'*r̃)[1:rnk] | |
mul!(view(p.delbeta, 1:rnk), view(p.qr.Q, :, 1:rnk)', r̃) | |
ldiv!(R, view(p.delbeta, 1:rnk)) |
?
Co-authored-by: David Widmann <[email protected]>
Avoid erroring out for low rank design matrices when dropcollinear=false.
Avoid unnecessary triangular solves. Avoid indexing in the Q.
Avoid slicing R matrix in a way that triggers a minimum norm solution.
Remove unnecessary temporaries.
Fixes #558
I removed the white space so please review with https://github.com/JuliaStats/GLM.jl/pull/559/files?w=1